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ABSTRACT 



This thesis deals with improving the miss distance of a missile, with imaging seeker(s), 
by utilizing dynamic image processing. In an encounter with a missile, a target tries to 
avoid the missile by performing an evasive maneuver when the missile is at a relative dis- 
tance which maximizes the miss distance. Dynamic image processing permits us to identify 
the evasive maneuver of the target by estimating its acceleration in magnitude and direc- 
tion. This thesis studies methods of utilizing this additional information about the target’s 
behavior in order to improve the missile’s performance. First the proportional navigation 
guidance law is explored in order to verify its advantages and weaknesses. Then, methods 
of obtaining the time dependent 3-D movement of a target from its image plane feature 
point correspondences are derived. The 3-D components of the target’s acceleration are ob- 
tained by using a Kalman filter. Missiles with two cameras, one camera and one seeker (ra- 
dar or IR), and only one camera are considered. Methods to get stereo vision by using the 
one camera plus one seeker setup and the single camera setup are proposed. Advanced 
guidance laws, namely advanced proportional navigation and optimal guidance are de- 
rived, for a 3-D environment. A three dimensional simulation program is developed using 
classical proportional navigation, advanced proportional navigation, and optimal guidance. 
The engagement is simulated using state variable design and the performance of the guid- 
ance laws is compared. 
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I. INTRODUCTION 



The U.S, as a result of the highly effective kamikaze attacks during World War II on 
U.S vessels, initiated the development of the first tactical missile (Lark guided missile). 
Since that time, proportional navigation guidance has been used in virtually all the world’s 
endoatmospheric tactical radar, infrared(IR), and television(TV) guided missiles. 
Proportional guidance works well not only for predictable targets, but also for highly 
responsive ones (i.e. targets executing evasive maneuvers). The proportional navigation 
guidance technology currently in use appears to be adequate, if the effective time constant 
of the guidance system is short in comparison with the flight time and, if the missile has 
considerable acceleration advantage over the target. The popularity of this interceptor 
guidance law is the result of its simplicity of implementation, and effectiveness. Although 
proportional navigation was apparently known by the Germans during World War n, no 
applications of it were reported. In the U.S, this guidance law was studied under the' 
auspicious of the U.S Navy. Proportional navigation was originally conceived from 
physical reasoning. The mathematical derivation of the “optimatility” of proportional 
navigation came more than 20 years later. 

This research develops a three dimensional missile/target simulation using three 
techniques of interceptor guidance, namely classical proportional navigation, augmented 
proportional navigation and optimal guidance. The primary research goal is to improve the 
miss distance of a missile with imaging seeker(s) by utilizing dynamic image processing. 
The existent dynamic image processing algorithms can be used to estimate motion 
parameters of the target. This additional information, about the target behavior, will be 
included in the proportional navigation homing loop in order to increase the missile 
percentage of kill by improving the final miss distance. Information about the target motion 
is especially important in the final phase of the engagement, given that an evasive 
maneuver performed by the target creates appreciable miss distance that may preclude a 
target kill. In an encounter with a missile, a target tries to avoid the missile by performing 
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a evasive maneuver when the missile is at a relative distance that maximizes the miss 
distance. A simulation of the adjoint model of the linearized homing loop permits us to 
obtain miss distance projections as a function of flight time or, if preferable, as a function 
of the time to go. The target can induce the most miss distance by executing an evasive 
maneuver at a short time to go. More precisely, the optimal evasion from the target “point 
of view” would be a series of maneuvers at the times of flight that, by superposition, 
produce the most miss distance. Estimating the target maneuver and incorporating this 
information into the guidance control input is perhaps the difference between success and 
failure. 

Chapter II introduces the idea of proportional navigation and how the actual guidance 
law is developed. Chapter HI deals with estimating the target motion parameters by using 
two perspective views. In Chapter IV, the augmented proportional navigation and optimal 
guidance will be derived. Also in this chapter, a tridimensional missile/target simulation is 
developed using classical proportional navigation. Subsequently, the target’s estimated 
motion parameters will be incorporated in the tridimensional engagement by using 
augmented proportional guidance and optimal guidance. Chapter V consists of actual 
simulation results. The different control laws will be tested and compared by producing 
miss distance projections for different evasive maneuvers. Finally, conclusions and 
recommendations follow in Chapter VI. All computer simulations are developed using 
Matrix Laboratory(MATLAB). 
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II. FUNDAMENTALS OF TACTICAL MISSILE GUIDANCE 



A. GENERAL 

Proportional navigation guidance (PROPNAV) commands the missile to turn at a rate 
proportional to both the angular velocity of the line of sight (LOS) and the closing velocity. 
The constant of proportionality is a unitless designer chosen gain (usually in the range 3-5) 
known as the effective navigation ratio/constant. Mathematically, the guidance law can be 
stated as 

u m = NV c X, (Eq 2.1) 

where u m is the acceleration command which is perpendicular to the instantaneous LOS. N 
is the effective navigation constant V c is the closing velocity along the LOS, and X is the 
LOS angle (in rad). The overdot indicates the time derivative. 

If the navigation ratio is greater than 1, the missile will be turning faster than the LOS, 
and thus the missile will build up a lead angle with respect to the line of sight. For a constant 
velocity missile and target the generation of this lead angle can put the missile on a collision 
course with the target (zero angular velocity of the line of sight). If N = 1 then the missile 
is turning at the same rate as the LOS, or simply homing on the target. If N < 1, then the 
missile will be turning slower than the LOS, thus continually falling behind the target, 
making an intercept impossible. In order to completely understand the physics of 
proportional navigation guidance it is necessary to analyze pursuit and constant bearing 
guidance. 

B. PURSUIT GUIDANCE 

For pursuit guidance, the missile velocity vector is always directed toward the target 
as illustrated by Figure 2.1. The missile is then constantly heading along the line of sight 
from the missile to the target and its path describes a pursuit path. Given that the rate of turn 
of the missile is always equal to the rate of turn of the LOS, “pure” pursuit (without leading 
angle) paths are highly curved. This requires the missile to use significant acceleration. 
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Since the signal processing is limited to continuously locating the target and changing the 
missile flight path angle, the on-board avionics are relatively simple. As will be 
demonstrated later, this kind of classical guidance law is a special case of PROPNAV when 
the effective navigation ratio is equal to 1 . 




Figure 2.2 shows the geometry of the pursuit guidance law. V m and V t are respectively 
the missile and target velocities, Q m and 0, are respectively the missile and target flight 
path angles, a f is the difference between the LOS angle and the target flight path angle, and 
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r is the instantaneous separation between missile and target. Inertial and missile translating 
coordinate systems are also shown in the figure. 

The velocity of the target with respect to the missile is given by: 



v = V,- V m \ (Eq 2.2) 

v = re r + r0e e . (Eq 2.3) 




Writing the velocity of the target and missile in terms of the polar unit base vectors e T and 
e % , we get: 
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(Eq 2.4) 



V, = VfCoscq*?,. - KjSina^Q, 

Vm = V m e r . (Eq 2.5) 



Equating equations 2.2 and 2.3 and using equations 2.4 and 2.5, we obtain: 



r = V lC osa r V m , 


(Eq 2.6) 


= ~ v m^ na r 


(Eq 2.7) 


From Figure 2.2, we see that: 


0 m = a + 0 

m t t 


(Eq 2.8) 


Considering a non responsive target: 


•s 

II 

£ 

•o 


(Eq 2.9) 



The missile acceleration u is obtained by differentiating equation 2.5: 

“ = ^r + ^/e> (Eq 2.10) 

given that from analytical mechanics: 



de r 






dt Sd dt ' 



(Eq 2.11) 

Assuming constant speed (magnitude of the velocity vector), the acceleration command 
will be the normal component of the acceleration which will be designated u m : 



u m = V m Q m = V m X = V m*r 



where d f is a time function. 



(Eq 2.12) 



C. CONSTANT BEARING GUIDANCE 

The accelerations required by the pursuit guidance law can be reduced by aiming the 
missile ahead of the target by using a lead angle. In this case the missile traverses a straight 
line to a collision with a constant speed non maneuvering target, as shown in Figure 2.3. 
The missile converges on the target by using a constant LOS angle ( X = constant). Since the 
rate of change of the LOS angle is zero throughout the flight, the lateral accelerations are 
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zero. If the target maneuvers evasively, by changing its velocity vector in direction and/or 
in magnitude, a new collision course must be computed and the missile flight path altered 
accordingly. The constant bearing geometry is shown in Figure 2.4. 

We wish to fmd the missile control input necessary to responde to prescribed target 
accelerations. The relative velocity is given by: 



V = V,-V m = re r + rXe Q . 


(Eq 2.13) 


The target and missile absolute velocities are: 


V t = V [ cosa [ e r -V [ sin (x t e e ; 


(Eq 2.14) 


V m = VUcosa e-V m sma e ft . 

rn m m r m m ” 


(Eq 2.15) 


Subtracting equation 2.14 from equation 2.15 and equating 


the result to 2.13, we find that: 


r = V t cos a-V m cos a m ; 


(Eq 2.16) 


rX = -K, since, + ^sina . 

t t m m 


(Eq 2.17) 


The requirements for a constant bearing guidance are: 


X - 0; 


(Eq 2.18) 


r <0. 


(Eq 2.19) 


Using equations 2.17 and 2.18, we get: 



sin a 

sina m = -rr-V,, (Eq2.20) 

v m 



From equations 2.16 and 2.19: 

cos a, 

cos a m > ~u — V t ' (Eq 2.21) 

v m 
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Figure 2.3 Constant Bearing Guidance Trajectory 



We may use equation 2.20 to obtain a expression for coscc f by squaring the equation and 
using the fundamental trigonometric entity. Doing this, we get: 



cosa, 

-V = 

V„ ' 



m 



v) 



m 



1 1 + (cosaj 2 



(Eq 2.22) 
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Substituting this expression into equation 2.21, we obtain: 



cos a > 

m 



(— 



This expression is satisfied if: 



1 

2 



+ (cos a ) 2 

x m' 



(Eq 2.23) 
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(Eq 2.24) 



V) 

-±-\<Q=>V m >V r 

* m 

Thus, for this guidance law to be effective the missile must have speed advantage relative 
to the target. 

The missile and target accelerations can be computed from Figure 2.5, which expands 
the acceleration in terms of tangential and normal components. The velocity vector of a 
body (missile or target) is described by: 

V = Vx, (Eq 2.25) 

where t represents the tangential unit vector to the trajectories represented by the dashed 
lines in the figure. The figure shows the response of the missile to a evasive target. We are 
interested in computing the relationship between the missile control input and the target 
maneuver. 




Figure 2.5 Constant Bearing Normal Acceleration 
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The acceleration of the body is given by: 



JV w ■ dx 

u = a -l = <L(Vx) = Vi + V— . 
dt dt dt 



(Eq 2.26) 



For small values of A\j/ it approaches the magnitude of Ax and the direction of Ax becomes 

dx 

perpendicular to the direction of x . It follows that the derivative is of magnitude 1 and 

d\y 



perpendicular to x . Then this derivative is the unit normal vector r'l . The time derivative 
is found by using the chain rule, as follows: 



dx 

dt 



dx dx d\\i 

— = = mjr = rid. 

d't dyd't 



(Eq 2.27) 

Assuming constant speed, the missile and target accelerations are always in the direction of 
the respective unit normal components and can be written as: 

(Eq 2.28) 

u, = n,V t O t (Eq 2.29) 






since. 



d = a m + X, 

m m 



(Eq 2.30) 



and given that X = 0, 



0 m = a m , 

m m ’ 



ii = V m a . 

m m m 



(Eq 2.31) 

(Eq 2.32) 

The variable, u m is the missile acceleration magnitude. Differentiating equation 2.20, we 
find d m and then u m in terms of the target evasive acceleration: 



d m {t) = 



cos a { ( t ) 



y m ^sa m (t) ' 

where the time factor is included. Additionally, 



a, (0 . 



(Eq 2.33) 
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. u ,(0 

a,(0 = 9,(0 = — 77 - 

v t 



(Eq 2.34) 



where u,(0 is the target acceleration magnitude. Hence: 

/ cosa f (;) \ 



u m( 0 = 



1 cos a (0 



u ,(0 • 



(Eq 2.35) 



From this last equation and equation 2.20, we conclude that the LOS will maintain its 
direction in space, keeping the missile on a collision course with the target provided that 
the missile’s and target’s kinematics normal to the LOS behave likewise. Additionally, 
from equation 2.21, the closing velocity (component of the relative velocity along the LOS) 
must be positive. Constant bearing guidance requires the knowledge of the heading and 
velocity of the target, the line of sight, and the velocity of the missile, which dictates a more 
complex signal processing system than for pursuit guidance. 



D. PROPORTIONAL NAVIGATION GUIDANCE 

1. In Search Of The Proportional Navigation Concept 

Pursuit guidance tries to continuously point the missile to the target, resulting a 
highly curved path and very large accelerations. The guidance law is only interested in the 
present position of the target; lacking information about the target kinematics. This lack of 
information precludes the missile from building a lead angle, resulting in a somewhat 
ineffective guidance law. Constant bearing guidance points the missile to the future 
position of the target, resulting in a straight line collision path with a non maneuvering 
target. Before pointing the missile, the guidance system needs to know the heading and 
velocity of the target to compute the target’s future position. So, this method is not 
practical, especially when dealing with targets with evasive capabilities. 

The advantage of proportional navigation is that it provides a practical method of 
approximating a constant bearing course to a maneuvering target. PROPNAV tries to 
emulate the constant bearing guidance command by using LOS rate information from an 
on - board electromagnetic or electro - optic device. 
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A missile using constant bearing guidance only needs a control input when it is 
necessary to change its heading at the beginning of the flight and afterwards if the target 
maneuvers. The form of this command signal was derived and is repeated here for 
convenience: 

n m = V mV- (Eq 2.36) 

From the proportional navigation geometry in Figure 2.6: 

e„ = a m + X. (Eq 2.37) 

Taking its derivative: 

= «„ + *•• (Eq 2.38) 

The acceleration of the missile (assuming constant speed) is: 

= V m 6 m » = V m (a m + X)h„ (Eq 2.39) 

Our goal is to emulate equation 2.36 by using a linear transformation between the LOS rate 
k and the missile’s angle rate a m . Set, for example: 

j.= J-T<V (Eq 2.40) 

Equation 2.39 becomes: 

1 N 

u m ~ ^m(« m + T 7 — rcOft = V m — — r dji. (Eq2.41) 

m m v m m _ 1 m' m vf _ i m 

By letting N be large this equation approaches equation 2.36 for a constant bearing path 
(collision course). We are interested in the relationship between the LOS rate and the flight 
path angle rate. Using equations 2 38 and 2.40, we find that: 

Q m = a m + k = (N- l)k + X = NX- (Eq 2.42) 

u m = NV m Xn. (Eq 2.43) 

Therefore, PROPNAV is a practical guidance law that emulates constant bearing guidance 
by issuing control commands that are proportional to the LOS rate. 
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Pursuit guidance is a particular case of proportional navigation when N =1 
(compare equations 2.12 and 2.43). As we have seen, constant bearing guidance is obtained 
by letting N be large (theoretically infinity). However, large gains in the amplifiers also 
cause large amplifications of noise; therefore N is usually restricted to less than 6. 
proportional navigation paths are less curved than pursuit paths, but more curved than 
constant bearing collisions. PROPNAV anticipates the future position of the target without 
actually computing it. Due to this property this guidance law presents a higher degree of 
responsiveness than other guidance laws. 

2. Proportional Navigation And Zero Effort Miss 

In Figure 2.6 the missile, with velocity magnitudeV m , is heading at an angle of 

L + HE with respect to the line of sight. The angle L is known as the missile lead angle and 
is the theoretically correct angle for the missile to be on a collision triangle, with the target. 
If the missile is launched in a collision triangle with a non evasive target, no further 
accelerations commands will be required to hit the target. The angle HE is known as the 
heading error, and represents the initial deviation of the missile from the collision triangle. 

In practice, the missile is usually not launched exactly in a collision triangle, since 
the expected intercept point is not known precisely. The location of the intercept point can 
only be approximated, because we do not know in advance what the target will do in the 
future. In fact, that is why a guidance system is required. The point of closest approach of 
the missile and target is known as the miss distance. Guidance system lags or subsystem 
dynamics will cause miss distance. The simplest proportional navigation homing loop is 
shown in Figure 2.7 where we have linearized the missile/target engagement by using the 
small angles approximation (i.e. we assume that the flight-path angles and the line of sight 
angle are small in order to linearize the engagement geometry. Then, the cosine functions 
are approximated by 1 and the sine and tangent functions by their arguments). In a 
linearized analysis, the range equation is approximated by the following time varying 
relationship: 
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(Eq 2.44) 



r = V c {t F -t) = V c t g0 
where V c is the closing velocity, t F is the total flight time, and t (time to go) is the time 
until the end of the flight. As shown in [Ref. 1] the miss distance will always be zero in a 
zero-lag proportional navigation homing loop. The PROPNAV guidance law used in the 
homing loop of Figure 2.7, and also the most used in the literature, is not the one derived 
in equation 2.43, but the following one: 

u m = NV c Xn x (Eq 2.45) 

where n x is the unit vector normal to the LOS. Then, the control input is issued 

perpendicular to the instantaneous LOS. It can be easily demonstrated that this last 
expression maintains the proportionality between the missile flight path angle and the 
angular LOS rate. In [Ref. 2], Guelman contrasted “pure” PROPNAV (described by 
equation 2.43, wherein command accelerations are normal to the missile velocity vector) 
and “true” PROPNAV (described by equation 2.45, wherein command accelerations are 
normal to the line of sight). He concluded that the later law would result in intercept only 
if the initial conditions were within a well-defined subset of the parameter space. In the 
homing loop of Figure 2.7, the seeker provides the LOS rate by taking the derivative of the 
geometric LOS angle. The noise filter processes the noisy LOS rate measurements to 
provide an estimate of the LOS rate. The guidance command is generated using the “true” 
PROPNAV guidance law. The guidance system must cause the missile to maneuver, by 
using moving control surfaces. The seeker and the guidance system dynamics are described 
by differential equations. 

The presence of delays in the homing loop creates miss distance. In the presence 
of guidance system dynamics, the heading error ( HE) and target maneuver (target evasive 
acceleration) are the two sources of miss distance. The PROPNAV guidance law can be 
expressed in terms of the zero effort miss. The zero effort miss is not only useful in 
explaining PROPNAV but is also useful in deriving more advanced missile control laws. 
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Figure 2.6 Proportional Navigation Two Dimensional Engagement 



The zero effort miss is the distance the missile would miss the target if the target 
continued along its present course and the missile made no further corrective maneuvers. 
Using Figure 2.6: 

ZEM X = r x + v x t g0 , (Eq 2.46) 

ZEM y = r y + v y t g0 ; (Eq 2.47) 

where ZEM represents the zero effort miss, r is the missile/target relative distance, and v is 
the missile/target relative speed. The subscript (x or y) represent the projection of the 
respective quantity over that coordinate axis. 
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Figure 2.7 Zero - lag Proportional Navigation Homing Loop 
(Linearized Engagement) 



The ZEM perpendicular to the LOS, is given by: 

ZEM pL0S = - ZEM X sin X + ZEM y cos X. 
Using equations 2.46 through 2.48, we obtain: 



t ( r v — r v ) 

7UAA _ X V y y X' 

ZEMpios 



The LOS is: 




taking its derivative, we obtain: 



X = 



r v — r v 

XV V X 



(Eq 2.48) 



(Eq 2.49) 



(Eq 2.50) 



(Eq 2.51) 
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Comparing equations 2.49 and 2.51, the LOS rate may be expressed in terms of the 
component of the zero effort miss normal to the LOS: 



X = 



zem plos zem plos 



rt 



go 



V c‘lo 



(Eq 2.52) 



where r = V c t ■ Then the PROPNAV guidance command magnitude can be expressed in 
terms of the ZEM perpendicular to the LOS: 

N ZEM . 



u m = 

m 



PLOS 



(Eq 2.53) 



go 



Thus, we conclude that the PROPNAV acceleration command that is perpendicular to the 
LOS is not only proportional to the LOS rate and closing velocity but is also proportional 
to the zero effort miss and inversely proportional to the square of time to go. The efficiency 
of PROPNAV guidance is a direct consequence of this dynamic property. This is, of course, 
a very powerful concept. 
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III. DYNAMIC IMAGE PROCESSING 



A. GENERAL 

A missile that uses a TV camera and a seeker (radar or IR), or instead, two TV cameras 
is considered. A setup with only one TV camera is also studied. The seeker and the camera, 
or the two cameras, can be located on the missile’s nose separated by a transversal distance 
d. The seeker plus the single camera setup, permits the missile to emulate the stereo vision 
of the two cameras setup. It has the additional advantage of tracking the target at the early 
stages of the engagement using solely the seeker’s LOS angle information. This system 
permits us to compute the 3-D target motion by using a two perspective views motion 
algorithm and the target’s spatial direction and range provided by the seeker. The two 
cameras setup permits us to use image plane locations in two views, corresponding to a 

single object point at times t, and t 2 , to determine the 3-D object (target) locations X 0 (z^) 

and X 0 (r 2 ) . The one camera setup also permits us to determine the motion of the target, as 

a function of time. This is done by using a two perspective views motion algorithm and 
guessing the target’s physical dimensions to estimate its absolute depth. In this way, we 
emulate binocular vision. The estimated 3-D motion of the target and the image sampling 
time permit us to estimate the target velocity and acceleration components in a preselected 
3-D rectangular coordinate system. The acceleration information can subsequently be 
injected into the control algorithms, which will be developed in the next chapter, to improve 
the miss distance. 

1. Scene (3-D) - Image (2-D) Geometric Considerations 

Mathematically, we can express the transformation of object point locations (3- 
D) to image plane locations (2-D) by the following generally noninvertible geometric 
transformation: 

Xi(t) =g(X 0 (/),...). (Eq 3.1) 
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The modeling of the imaging process, described by the above equation, relates object points 

X 0 (t) in the 3-D scene to image points X,(i) in the image plane. The function g depends 
on the imaging geometry, lens model, and coordinate system choices. 

The imaging model is derived by considering the pinhole camera model shown in 

Figure 3. 1. The point X c lies over the camera's optical axis at a distance /from the image 
plane. The figure shows two distinct coordinate systems, an image plane coordinate system 
and a global coordinate system. Our first goal is: assuming the simplified camera model 
shown in Figure 3.1, derive the transformation described by equation 3.1 where the object 

point X 0 (t) can be measured from either coordinate system. 
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The object/image relationship defined in equation 3.1 is defined by a transformation 
matrix. Independent of the camera model, this transformation matrix is the product of two 
matrices. The first matrix describes the object - image coordinates transformation, and is 
derived by assuming that the 3-D object coordinates are measured relatively to the image 
plane. However, if the object coordinates are measured relatively to the global coordinate 
system, a second transformation matrix relating the two coordinate systems have to be 
defined. This matrix is the composite of the relative rotation and translation between the 
coordinate systems. It describes the coordinates transformation between the two coordinate 
systems. 

To identify the transformation defmed by equation 3.1, the two matrices are 
derived for the simplified camera model of Figure 3.1. Monocular vision (only one camera) 
is incapable of determining absolute depth. However, any imaged point is constrained to 

correspond to an object point located anywhere on the 3-D line segment containing Xj (r) ,• 



X- and X Q (r) . 

Assuming that the coordinate systems for both object and image points, are 
coincident and centered in the image plane, the above colinear points are related by: 



*(X/(0 -X c ) = (X c -X 0 (t)) . 



Expanding this equation yields: 



^[ 0 ^ z 3 r "^ 00 ] r} 




(Eq 3.2) 



(Eq 3.3) 



where the superscript T is the transpose operator. The time index t has been dropped to 
simplify the notation. Equation 3.3 yields: 



, ( x o -f) x o , 


(Eq 3.4) 


/ f l ' 


y 0 fy Q 


(Eq 3.5) 


yi ~ k " (f-x 0 )' 
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7 . 

“ l 



(Eq 3.6) 



“o _ J~o 

1 " 

The minus sign in the second expression of the two last equations, stands for the image 
inversion originated by the back - projection model of Figure 3.1. The matrix representation 
of the nonlinear equations 3.5 and 3.6 is: 



Xi 



MX 0 = 



0 / 00 
0 0/0 
L-l 0 0/J 




(Eq 3.7) 



The last equation uses homogeneous coordinates (a technique also used to develop 
computer graphics), for image and object points. The homogeneous coordinates are defined 
by multiplying the physical coordinates by an arbitrary constant c and including the 
constant as an additional element of the vector: 




n r 



cy ( - cz ( - c , 



and 



(Eq 3.8) 



r i T 

x o= x 0 y 0 z o ] - ■ 



(Eq 3.9) 

Note that in equation 3.9 the arbitrary constant is equal to 1. The object - image point 
transformation is defined by equation 3.7, where it is implicitly assumed that x i = 0. 
Rewriting equations 3.5 and 3.6, we conclude: 



y_o 

y-t 



z i 



x o~f 
-f • 



(Eq 3.10) 



Fixing the image plane coordinates y i and z ( - the above equation describes the 3-D line over 

which the 3 -D object is located. Therefore, while the transformation in equation 3.7 is not 
invertible, choice of a specific image point constrains corresponding object points to lie 
along a 3-D ray (shown in Figure 3.1). 

If the object points are measured relatively to the global coordinate system, the 
matrix relating the two coordinate systems has to be computed. This transformation matrix 
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is the product of a succession of matrices. Individually, each of these matrices defines a 
rotation or translation of the image plane coordinate system relative to the global coordinate 
system. The succession of transformations may be of the form: 

Xo = ( T2 R2 Rl Tl )X*o> (Eq 3.11) 



where X 0 and xf define the homogeneous object coordinates in the image plane and global 
coordinate systems, respectively. Here the first transformation is the translation Tl 
followed by the rotation Rl, etc. The composite of the above transformation may be defined 
by: 



H g ^i= T2R2RITI- (Eq 3.12) 

Then, the general relationship between object points measured relatively to any user 
selected coordinate system and the image plane points is: 



X,- = (Eq 3.13) 

For the simple case of only a translation as shown in Figure 3.1, we see that in object 
coordinates: 



*0 = l-\d x d y d : 



(Eq 3.14) 



Homogeneous coordinates enable us to represent the last relationship using a translation 
matrix: 



X 0 = H g . 



;Xi = 



,10 0 
0 1 0 -d y 
0 0 1 -d : 




000 1 



(Eq 3.15) 



2. Stereo Vision (2 Cameras) 

Monocular vision disables depth perception. In fact, due to the impossibility of 
inverting the 3 x 4 matrix Q = MH g • obtained in the last section, we are constrained to 
the determination of image points from object points. However, we are interested in 
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determining the 3-D locations (measured relatively to a global coordinate system). One 
approach to solve this problem is to use more than one camera. One of our proposals, is to 
use two cameras in the missile’s nose separated by a distance d emulating, in some way, 
the human visual system. 

Initially, we assume the simplified two dimensional diagram of the stereo vision 
in Figure 3.2. The scene consists of a 2-D surface. As shown in the figure, a point on this 
surface is projected onto the two image planes (IP1 and IP2). In general the two centers of 
projection differ in length (/j and / 2 ). It is assumed that the user selected global coordinate 
system, to measure the object coordinates, is coincident and centered in the image plane 
IP1. The coordinate x i is shown in the figure, the coordinate y i is perpendicular to x t and 
in the plane of the page. The object point may be determined using the two projected points, 
one in each camera. 

The relationship between the homogeneous coordinates of the object point, 
measured relatively to the image plane EP1, and the homogeneous coordinates of the 
corresponding image point is: 



cx n 

c 



1 

o 

O, 




*0 


0 ~T 1 






L a _ 




_1_ 



(Eq 3.16) 



The object coordinates measured from IP2 are related to the object coordinates measured 
from IP1 by the following relationship: 




or, in homogeneous coordinates: 




1 0 d 
0 1 0 
0 0 1 



1 

A 0 



(Eq 3.17) 



(Eq 3.18) 
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Figure 3.2 Two Dimensional Stereo Vision 



Then image points in IP2, denoted x i2 , may be related to object points measured relatively 
to the global coordinate system by: 



r 




"l o o 


1 0 d\ 


x 0 


cx i2 


= 


0 1 


0 1 0 


y 0 


c _ 




L h \ 


0 0 1_ 


X 



(Eq 3.19) 



Using equations 3.16 and 3.19 the following relationships in object coordinates are 
obtained: 
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*«i C/j -y 0 ) 



(Eq 3.20) 



* 






o 



xnVi-yo) 
f 2 



- d. 



(Eq 3.21) 



Equating these two equation the object’s depth v 0 may be found: 

_ f\fi ( x n ~ x a ~ d) 
y ° f\ x n ~fi x i 



(Eq 3.22) 



The object point x 0 may be found, from either equation 3.20 or equation 3.2 1 . Hence, using 
two image planes permits us to determine the object point depth y 0 from its corresponding 

image points. This was proved for a 2-D surface. Next we are going to see how to do it in 
a 3-D environment. 

Equation 3.13, defines the relationship between the scene three dimensional 
points and the correspondent two dimensional image points by using a matrix of 
transformation: 



Q = MH g _^ i . (Eq 3.23) 

H g • is the coordinate systems transformation matrix which depends on the rotations and / 
or translations of the image plane coordinate system relative to the user selected global 
coordinate system. M is the object - image transformation matrix, which is a function of 
the imaging geometry and lens model. The Q matrix is a non - invertible 3x4 matrix 
(assuming homogeneous coordinates) and may be generically represented by: 



Q = 



^11 ^12 ^13 Q 14 
<?21 <?22 *?23 24 • 

531 ^32 ^33 34 



(Eq 3.24) 



The missile must have sufficient processing capability to fmd this matrix in real time. The 
object - image transformation matrix M is generally invariant (however zooming the scene. 



for example, changes its value). The coordinate systems transformation matrix H g ( must 
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be dynamically updated as the missile/target engagement proceeds. The selected global 
coordinate system for missile guidance simulation purposes is a ground coordinate system 
which will be presented in the next chapter. For full dimension (3 -D) stereo vision, the two 
cameras arrangement may be described using homogeneous coordinates as: 

X ic = Q C X 0 , (Eq 3.25) 

where the index c = 1, 2 refers to the cameras. Since each of these two matrix equations 
(one for each sensor) represents two equations in physical coordinates, we obtain four 
equations and three unknowns by using the two cameras stereo arrangement. The matrix 



(Eq 3.26) 



(Eq 3.27) 

The index of each matrix element is composed by three numbers, the first is the camera 
number and the next two represent the element position into the matrix. Each of these two 
matrix equations generates two equations in physical coordinates. To find these equations, 
the arbitrary constants (Cj and c 2 ) have to be calculated. Then, each of the constants is 
substituted into the two remaining matrix equations. Finally, regrouping terms as 

coefficients of x g , y g 0 and z g , a set of four equations is obtained. Performing this procedure 
to obtain the first physical equation from the matrix equation 3.26, we get: 

C 1 = <7 131 4 + <?132 y 8 o + *?133 z o + *7134’ 3 - 28 ) 



equation 3.25 can be explicitly written for each sensor as: 



and. 
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*7m *7ii2 *7 113 *7 n4 


c l z «'l 




<7 121 *7 122 *7i23 *7 124 


_ C 1 




_*7l31 ^7 132 *7 133 *7 134 


C 2- V i2 
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C 2 Z «2 


ZZ 


*7221 *7222 *7223 *7224 


L C 2 _ 




*7231 *7232 *7233 *7234J 
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substituting this expression into Cjy ( j, and regrouping terms, we obtain: 



(Eq 3.29) 



((<7in- <7131^1) 4+ (<7i2i-<?i32. v a). v o + (4n3 ~ ^133^il) z o ~ < 7i343',i “^iu) ' 



The set of four equations and three unknowns is written compactly in matrix notation as: 



PXl = F 

or 



(Eq 3.30) 



*7111 — *7131^11 7 121 — 7 132^ il *7 1 13 — 7 133^ il 
7121 -7 131 2,1 7 122 — 7 132“ 1 1 7l23 - 7i33 z il 
7211 — 7231^12 7212 “ 7232^2 7213 — 7233> 7 i2 
7221 “ 7231 z i 2 7222 " 7232 z i2 7223 “ ^223 Z i2_ 



X 8 

A o 



7 134 X il 7ll4 

7134 Z il" 7124 (Eq3 31) 
7234 - v i2“ 72141 
_7234 z i2 ” 7224! 



Equation 3.31 may be solved using least square techniques by forming the pseudoinverse 



of P denoted Ft . Hence: 



PX 8 0 = F^P T PX 8 0 = P T F=>X 8 0 = ( P T P)~ l P T F=>X 8 0 = PfF. (Eq 3.32) 

This equation yields the mean square estimate for the object point X 8 0 . Alternatively, X 8 0 
may be found by using three of the four equations, assuming that the three equations are 
linearly independent. 

In this exposition, we have assumed that the necessary image plane point 
correspondences have been determined. It is important to say that this is the most difficult 
problem in the development of a stereo vision algorithm. Techniques to solve this problem 
are presented in [Ref. 3]. [Ref. 4] presents an algorithm to match stereo images. 
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B. ESTIMATING 3-D MOTION PARAMETERS OF A RIGID BODY FROM 
TWO CONSECUTIVE IMAGE FRAMES 



1. General 

In section A of this chapter, we have shown that using two cameras, we are able 
to find the object (target) 3-D position relative to a user selected global coordinate system. 
As intermediate steps, it is necessary to establish feature correspondences between selected 
points in the two stereo images (static stereo). It is also necessary to establish feature 
correspondences for all pairs of consecutive image frames in each camera’s image 
sequence. The targets that we are interested in are mainly airplanes. Therefore, we may use 
as points for feature correspondences the tips of the wings, nose, stabilizers and rudder. 

Estimation of the 3-D target acceleration components may be divided in three 
steps. In the first step, the target estimated points at t l and t 2 are used to estimate its 3-D 
velocity components. In the second step, the 3-D velocity components of the target are 
computed using the target’s estimated points at t 2 and / 3 . Finally, the third step estimates 

the target’s 3-D acceleration components by identifying the time change of the target’s 
velocity for each of the three velocity components. Alternatively, we may use Kalman 
filtering theory to estimate the 3-D target acceleration components. 

A different approach for estimating the 3-D, time dependent, target motion is now 
presented. The formal definition of a rigid 3-D object is one for which the 3-D distances 
between any pair of points on the object do not change with time that is, for all pair of points 
on the 3-D object: 

| X m - X J 2 = c mn , Vf, V/n, n ; (Eq 3.33) 

where c mn are constants. The assumption of an rigid, or nondeformable target is reasonable 

and creates additional constraints for motion estimation. Rigidity constrains the motion of 
individual object points to be strongly coupled, although the need for point correspondence 
maintains. Then, the 3-D translation of the target may be determined by estimating the 
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translation parameters of a single point object. The basis of estimating the target’s motion 
using this approach, is that the 3-D motion of a rigid body can be described by a 3-D 
translation vector and three rotation angles chosen with respect to a user selected coordinate 
system. Then, six parameters completely defme the target’s motion. Formulating the 
rotations using three rotation matrices (/?g , /? a and R^), the target’s motion is described 



by: 



X 0 (t 2 ) = RX 0 (t j) + t. 



where R is the overall rotation matrix: 



(Eq 3.34) 



R = R e R a R p , (Eq 3.35) 

and T is the translation matrix. Equation 3.34 may be represented in homogeneous 
coordinates as: 





0 0 0 | lj 



(Eq 3.36) 



2. Monocular Motion Estimation Using Two Perspective Views 

Our goal is to compute dynamically the rotation (/?) and translation (7) matrices 
from point (or feature) correspondences between two perspective views. We can divide the 
process of estimating the three - dimensional motion of the target from image sequences in 
three steps. The first step is to establish feature correspondences between two consecutive 
image frames. Correspondences between features may be established through matching or 
inter - frame tracking. [Ref. 4] develops a two - view/stereo matcher that computes 
displacement fields from two images. The second step is to estimate the motion parameters 
(R and T matrices). The third step is to estimate the 3-D motion of the target using equation 
3.34. There are a number of papers in the digital image processing and computer vision 
literature dealing with estimation of the motion parameters of the 3-D motion of a rigid 



30 



body from two consecutive image frames. Weng, Huang and Ahaja [Ref. 5], propose an 
algorithm that given 8 point correspondences solves for a intermediate matrix called the 
essential parameter matrix (E). Then the Rotation matrix (/?) and the translational direction 

T 

(the unit vector ) are obtained from E. The magnitude of the translational vector ( 1 T ) 

and the absolute depths of the object points ( x ok and x' ok where x ok is the absolute depth 
of the object’s k feature point and x' ok is the absolute depth of the object’s k feature point 
after being rotated and translated) cannot be determined by monocular vision. The 

x ok x 'ok 

algorithm also solves for the relative depths (— and The algorithm is unable to 

estimate the target’s position because it is not possible to calculate the absolute depth. This 
agrees with intuition, due to the lack of invertibility of the 3-D to 2-D image transformation. 
To overcome this problem, we propose to estimate the absolute depth by correlating the 
dimensions of the target over the image plane with the guessed physical dimensions of the 
target. A relatively easy trigonometric approach permits us to estimate the target’s depth 
given the target’s physical dimensions. Another approach for emulating stereo vision is to 
use range and directional spatial information from the seeker. This information is combined 
with the monocular vision equations in order to estimate the target’s 3-D motion. 

In conclusion, monocular vision may be applied to estimate the motion of the 
target if additional information about the target is available (or guessed). The target’s 
acceleration components may now be estimated and injected into the missile’s control 
algorithm. 

3. Stereo Motion Estimation Using Two Perspective Views 

As we have seen, two or more spatially distributed sensors enables determination 
of the 3-D target motion. Here, our goal is to determine, not only the target’s motion but 
also the rotation and translation matrices that describe the motion. Given two consecutive 
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time samples or “frames” in each sensor, as well as image pomt correspondences, we may 
write: 



x }(0 = Q\^iK( 0. 

x](t) = Qg_+jX 0 (t ) , 



(Eq 3.37) 
(Eq 3.38) 



*J(' + 7,) = Ql-n*o(‘ + T M ) = Qj J 

0 0 0 



T I 

X 0 (t ) , (Eq3.39) 



1 



+ = Ql^oi' + T,) = Q]^\ 

- - - 

L 0 0 0 



I T 

I X 0 (t) . (Eq 3.40) 

| 

I 1 



The matrix Q l g _ > ,• is the, already known, 3-D to 2-D transformation matrix (the superscript 
refers to the sensor). T s is the time between two consecutive images. The system is 
represented in homogeneous coordinates. Assuming n point for point correspondences, a 
total of 8n equations in physical coordinates is obtained. However, the number of 



unknowns is 12 + 3n (9 elements of R , 3 elements of T and 3 elements for each X 0 ). This 
yields the constraint on the number of point for point correspondences: 

12 + 3/i ^ 8n. (Eq 3.41) 

Since n must be an integer, n'Z 3. Thus, three corresponding image points from two views 
in two frames are sufficient to determine both the motion parameters and the 3-D location 
of the object points. 

The three target acceleration components may be computed by taking the second 
derivative following the filtering of the target’s motion data. Alternatively, the target’s 
motion may be processed by Kalman filters to estimate its acceleration components. 
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This additional information about the target’s behavior may be used to improve 
the missile guidance towards the target. In order to effectively use this information, we have 
first of all to determine control laws that can use and produce better results if this 
information is available. This will be stressed in the next chapter. 



33 



IV. SIMULATION DEVELOPEMENT 



A. CONTROL ALGORITHMS DEVELOPMENT 

Chapter II introduced proportional navigation guidance. In this section we derive more 
advance guidance laws. Contrary to PROPNAV, these guidance laws use the estimated 
acceleration of the target as a additional input to the homing loop. As will be seen in 
Chapter V, the advanced guidance laws relax the interceptor acceleration requirements and, 
in general, yield smaller miss distances. 



1. Augmented Proportional Navigation 

Proportional navigation issues control commands that are proportional to the 
predicted zero effort miss normal to the line of sight ( ZEM pL0S ). That is, the missile 
guidance system tries to minimize the final miss distance between the target and the missile 
by issuing acceleration commands that are proportional to the miss distance, that would 
result if the missile made no further corrective acceleration and the target did not maneuver. 
Therefore, if the target maneuvers evasively it generates additional miss distance that is not 
accounted for in the PROPNAV guidance law. Augmented proportional navigation also 
issues guidance commands that are proportional to the predicted miss distance. However, 
for augmented PROPNAV, the miss distance is estimated by taking into account the 
maneuver of the target (target acceleration). The augmented PROPNAV target's 
acceleration dependent term will be calculated. This term is injected into the homing loop 
to enhance the guidance performance. In the following analysis, we follow the 
nomenclature of Chapter II and the geometry of Figure 2.6. 

The x component of the miss distance, for an evasive target, is computed as 
follows (the y component is computed similarly): 

(r,(0) = »,(0, (Eq 4.1) 

where t' represents time. Then, 



34 



h 

ZEM x {t ) = rj f = r x (t) + {v x (f)df. (Eq 4.2) 

r 

Where, ZEM x {t) is the x component of the zero effort miss, predicted at time t’ = t and 
r x (t) is the present missile - target relative distance along the x axis. But, 

<(0 = Jp (MO), (Eq 4.3) 

where ci x t (?) is the x component of the target acceleration. Then, 

Mr) r 

| dv x (f) = \a x ( {t")dt", (Eq 4.4) 

V, (0 t 

where t" is the variable of integration. Hence: 

V 

MO = MO + Jo? (/")<*". (Eq 4.5) 

Substituting this equation into equation 4.2, we get the expression for the predicted x 
component of the ZEM at time t: 

hr 

ZEM x (t) = r x \ tmtF = r x (t) + vjt) t g0 + (f) dt" dt" . (Eq 4.6) 

i i 



Where, t g0 = t F - 1 is the time to go. The y component of the ZEM (r) is obtained using 
the same reasoning, and is: 

hr 

ZEM y (t) = r y \ (mt ' = r y (t ) + v y (r) t g0 + JJaf (f") dt" dt " . (Eq4.7) 

f t 



r v 

The two interior integrals J 'a x ( (t") dt" and ja y t ( t " ) dt" are time functions; call: 

/ t 



and 



MO = ft(t")dt". 



(Eq 4.8) 
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yo = \a y t (ndf 



(Eq 4.9) 



Then, equations 4.6 and 4.7 may be expressed as: 



(Eq 4.10) 

Z£M,«) = r - M'> + >V (')<*„ + )df = r,(r) +v,«)< so + /i,(( <( ,) 



'f 



(Eq 4.1 1) 



ZEM y (t ) = /y = r v (r) + v v (r) r^ 0 + J/: v (f”) Jr" = r v (r) + v v (r)r g0 + /i v (r g0 ) ; 



where h x ( y o ) and h y ( t g0 ) are t g0 and maneuver dependent functions. The component of the 
ZEM that is perpendicular to the LOS, ZEM pL0S , is: 

ZEM PL0S (t ) = ZEM y (r) cos (Mr) ) -ZEM x (t) sin (MO ) . (Eq4.12) 



Substituting 4.10 and 4.11 into this equation and setting cos (Mr)) 



yo 

MO 



and 



sin (Mr) ) 



yo 

TUT 



we obtain: 



ZEMpiosi 1 ) - 



(Eq 4.13) 



r x (0 

(r y (r) +Vy (r)t go +h y (t go )) T ^- 



r (r) 

- (r (r) + v (r) t +h (t )) 

v x v/ x v/ go x v go /7 I* f A 



From equation 2.52: 



MO 



r x (t)v y (t) -r y {t)v x (t) 

r 2 (0 



(Eq 4.14) 



Then equation 4.13, which takes into account the target’s acceleration to estimate the miss 
distance, may be reduced to: 
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r x {t) r (t) 

ZEM PL0S (0 = t go r(t)\{t) + (h y (t g0 ) — -h x (t go ) . 



(Eq 4.15) 



The PROPNAV guidance command may be expressed in terms of the ZEM as (see derivation 
in Chapter II): 



"m (') 



W C (QZEM PLOS (0 
r <')'*. 



(Eq 4. 16) 



where ZEM pLOS {t) was then derived for a nonresponsive target. If the target maneuvers the 



zero effort miss is augmented by an additional term, on the right hand side of equation 4.15. 
Therefore, a perfectly plausible guidance law, in the presence of target maneuver, would be: 

(Eq 4.17) 



(') = NV C (0 X ( 0 + (hy ( , go ) cos ( X (,) ) - h x (t g0 ) sin ( X ( 0 ) ) • 



This guidance law is PROPNAV with an extra term that accounts for the maneuver of the 
target. The equation was derived for a nonlinearized geometry. The impossibility of knowing, 
a priori, the future target maneuver, precludes the calculation of h v (t on ) , h r (t on ) and t on . 
However, if the target desires to inflict the most miss distance it must maneuver at a small time 
to go. Also, considering the time constants associated with the target’s maneuver, we propose 
to approximate the time dependent target acceleration a t (t') by a constant target 
acceleration a t (t ) ) i.e. assume a,(t') = a t (t) . Then, the control law may be stated as: 

(Eq 4.18) 



m 



APN 



(0 = NVAt)\(t) + 



NV c (t) (ay(t)r 
r <')'«» 



go 



cos (X (0 ) - 



a x ,{t)r. 



go 



sin ()i(f)) 



or, 



i, = nvX + 

m APN c 



NVj 



c go 



2r 



(a y t cos(X) -afsin(X)) 



(Eq 4.19) 
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where the time factor was dropped. The extra term, present in augmented PROPNAV 
(expression in parenthesis on the right hand side of equation 4.19), is proportional to the 
component of the target acceleration normal to the LOS. 

Linearization of the nonlinear missile - target geometry is shown in [Ref. 1] to be 
an accurate approximation to the actual geometry. Then, assuming a linearized geometry, 
equation 4.19 may be further reduced to: 

u m = NV c X + ^a y - (Eq 4.20) 

where we have considered that the LOS angle is small. A zero - lag augmented proportional 
navigation homing loop assuming linearized geometry is shown in Figure 4.1. 
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Figure 4.1 Zero - lag Augmented Proportional Navigation Homing Loop 
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The additional target maneuver term required by the guidance law, appears as a 
feedforward term in the homing loop block diagram. 

2. Optimal Intercept Guidance 

The missile - target engagement scenario may be described in state space 
representation by the following linear system: 

J?(0 =/[.t(0,K(0.f] = Ax(0 +Bu(t); (Eq 4.21) 

.V is the n - dimensional state vector describing the relative movement between the missile 
and the target and also, the dynamics of the guidance system. The variable u is the m - 
dimensional missile’s control input vector. We seek to find a guidance law that is a function 
of the system states. There is an infinite number of possible guidance laws. Thus, it is 
necessary to state in mathematical terms what the guidance law should do. Certainly we 
wish to design a terminal controller that would bring certain components of x ( t F ) to zero, 

using “acceptable” levels of control. One way to do this is to minimize a performance index 
made up of a quadratic form in the control: 

t F t F 

J = ^L[x(t),u(t),t]dt = -fu 2 (t)dt; (Eq 4.22) 

o o 

subject to the terminal constraint: 

x t (t F ) =0, i = 0, 1, (Eq 4.23) 

and the constraints: 

x(t) =f[x(t),u(t),t] =Ax(t)+Bu(t), (Eq 4.24) 

x (0), given. (Eq4.25) 

In equation 4.23, p <.n. 

The miss distance will always be zero in a zero - lag PROPNAV navigation 
homing loop. Guidance system lags or subsystem dynamics will cause miss distance. 
Optimal guidance eliminates miss distance by canceling out the guidance system dynamics. 
In this way the optimal guidance law attempts to make the real world guidance system 
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appear to be a “perfect” (zero - lag) guidance system. To find the optimal control vector, 
u (t) , that brings the system from a initial state x (0) to a terminal state x (t F ) (where 

some of its components are zero), we can use the method of the Lagrange multipliers. Then 
the constraints (4.23) and (4.24) may be adjoined to the performance function (4.22) by 

using the multipliers e = (e, , £ p , 0 p+ , ..., 0 n ) T and A. = (X ( , X n ) T as follows: 

i ' f 

J = e T x(t F ) + - J (n 2 ( t ) + X T (t) [/Lv(/) + Bu ( t ) -x(t) ] } dt. (Eq 4.26) 
o 

The Hamiltonian is defined as follows: 

H[x(t),u(t),t) = L[x(t),u(t),t] +\ T f[x(t),u(t),t] . (Eq 4.27) 

Integrating the last term on the right hand side of equation 4.26 yields: 

(Eq 4.28) 
b 

J = e T x(t F ) -X T (t F )x(t F ) +X T (0)x(0) + J {H [jc(0, u(0,t] +X(t)x(t) } dt- 



Considering the variation in J due to variations in the control vector u ( 0 , we get: 

(Eq 4.29) 



dJ 



= [{e T -\ T )dx) tmtf + [X T dx] lm0 + l\(^ + X T (t))dx + ^ ( du dt 



dH 



In order to make the variations in J due to variations in u(t) independent from the 
variations in x ( t) produced by the variations in u ( t) , we choose the influence functions 
X (t) to cause the coefficients of dx to vanish: 



i r (0 = = 

ox dx ox 



Then: 



X(t) = -A T X(t), 



with boundary conditions: 



X(tf r) — £. 



(Eq 4.30) 

(Eq 4.31) 
(Eq 4.32) 
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Using these results, equation 4.29 becomes: 



dj = A r (0)t/.t(0) + ^dudt. (Eq 4.33) 

0 

Hence, X T (0) is the gradient of J with respect to variations in the initial conditions, while 
holding u (t) constant and satisfying the constraints of the problem. For an extremum, dJ 
must be zero for arbitrary du{t ) . This can only happen if: 




0<.t<,t F , 



(Eq 4.34) 



or: 



u T + X T B = 0. (Eq 4.35) 

Then, we may determine the control vector u(t) , as: 

u(t) = -B T X(t). (Eq 4.36) 

Substituting equation 4.36 into equation 4.24 and repeating equation 4.31, the following 
two - point boundary value problem is obtained: 



[A! 



4 -bb? 
o - a t ^ 



w 

L^j' 



The 2 n boundary conditions are: 

x(0), given. 



x i (^) — 0, 

Xj (tp) = 0 , 



i = l, 

i = p + 1, .... n. 



(Eq 4.37) 



(Eq 4.38) 
(Eq 4.39) 
(Eq 4.40) 



The n boundary conditions 4.39 and 4.40 may be replaced by the boundary condition 4.32, 
which may be rewritten as: 



\(t F ) = e.; i = 0 p , (Eq 4.41) 

\(r F ) = 0; / =p+l,...,n . (Eq 4.42) 
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The two - point boundary value problem 4.37, 4.38, 4.41 and 4.42 may be solved by the 

sweep method [Ref. 6], The sweep method seeks to find solutions of the form: 

MO = W(0JC(0 +T(0s; (Eq 4.43) 

z = U(t)x(t) + V(t)e, (Eq 4.44) 

where W (t) ,Y (t) , U (t) and V (t) are time dependent matrices, 

z r = [.t., ...,.v„] 1 and s r = [e,,...,e n ] = [X,,...,XJ . Therefore, we want 

1 P t~t F 1 P 1 P t m t F 

to find solutions for the influence functions X(t) that are function of the state vector x(t) 

and the final value of the influence functions, or equivalently, of the specified final states 

z. Since equations 4.43 and 4.44 must be valid at t = t F : 

W{t F ) = 0, (Eq 4.45) 

V(t F )=0, (Eq 4.46) 



Y(t F ) 



0 



pxp 



( n-p ) *p 



nxp 



U(t F ) 



Ipxp Opx {n-p) 

px n 



(Eq 4.47) 



(Eq 4.48) 



where I is the identity matrix and 0 is a zero matrix with the specified dimensions. 
Substituting equation 4.43 into 4.37 and treating s as a constant vector, we get: 



Wx + Wx + Ye = -A t ( Wx + Ye) . (Eq 4.49) 

Substituting jc from equation 4.37 into the last equation, and again using equation 4.43 to 
eliminate X, we obtain: 



(W+WA+ A T W - WBB t W ) x + ( A t Y +Y- WBB t Y ) e = 0. (Eq 4.50) 

This expression must be true for any x and e, so the coefficients of x and £ must vanish: 
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(Eq 4.51) 



W + WA + A T W - WBB r W = 0; W(t F ) = 0, 

and 

I p x p 

Y + A T Y-WBB T Y = 0; Y(t F ) = _ _ _ (Eq4.52) 

®n-pxp 

- 1 nx p 

Next, we differentiate equation 4.44 with respect to time, treating e and : as constant 
vectors: 

Ux + Ux + Ve = 0. (Eq 4.53) 

Substituting x from equation 4.37, and using equation 4.43 to eliminate X, we obtain: 

(0 + UA-UBB T W)x+ ( V-UBB T Y)e = 0. (Eq 4.54) 

This last expression must be true for any x and e, so the coefficients of .v and e must - , 
vanish: 

U+UA- UBB T W = 0, (Eq 4.55) 

V-UBB t Y = 0. (Eq 4.56) 

From equations 4.52 and 4.55 and the boundary conditions 4.47 and 4.48, we conclude that: 

U(t) = Y T (t). (Eq 4.57) 

Then equation 4.56 may be rewritten as: 

V = Y t BB t Y; V(t F ) = 0. (Eq 4.58) 

The Ricatti equations 4.51, 4.52 and 4.58 may be integrated backwards from the final 
conditions to yield W (/), Y (t) and V (t) . The equation 4.44 is solved for e to yield: 

e = [7(r)]- 1 [z-T r (r).v(r)]. (Eq4.59) 
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Our goal is, to find the influence functions A (/) in order to find the optimal 
control vector using equation 4.36. Equation 4.59 may be substituted into equation 4.43 to 
find A (t) : 



A (t) = (W-VV~ l Y r ).x(t) + YV~'z. 



However, the Ricatti equation 4.51 has as solution: 

W(t ) = 0. 

Hence, equations 4.52 and 4.58 become simply: 

Ipxp 

Y + A r Y = 0; Y (t F ) = - - - 

n 

v n-pxp 

- nxp 



and 



(Eq 4.60) 
(Eq 4.61) 



(Eq 4.62) 



h 

V(t) = -J {Y T BB T Y)dt. (Eq 4.63) 

t 

Combining equations 4.36, 4.60 and 4.61 yields: 

u(t ) = -B T X{t) = (-B T YV~ l ) [z-Y T x(t)] . (Eq 4.64) 

The final condition that we are interested on is z = [x , ...,x ] T = [0, ...,0]^. 

1 r » = » f 

Hence the expression 4.64 may be reduced to: 

u(t) = B r Y\r l Y r x(t), (Eq 4.65) 

where Y and V are computed from equations 4.62 and 4.63, respectively. 

Now that we have derived the terminal controller optimal feedback guidance law, 
we proceed to derive the continuous feedback law described by equation 4.65 for a single 
lag guidance system. The single - lag guidance system is mathematically described in the 
Laplace domain as: 
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a 



m 



U 



m 



1 

TTTs' 



(Eq 4.66) 



where a m is the missile’s acceleration, u m is the command acceleration, and T is the 

effective guidance system time constant. 

The relative motion between the target and missile is considered with the 
linearized (small angles approximation) intercept geometry shown in Figure 4.2. The 
assumption of small angles (flight path angles 0 , 0 f and LOS angle X) permits us to 

express the equations of motion in terms of state variables normal to the reference intercept 
course. The single - lag guidance model shown in Figure 4.3 integrates the missile - target 
relative motion of Figure 4.2 with the dynamics of the guidance system. The diagram of 
blocks of Figure 4.3 can be expressed in state space form as: 




(Eq 4.67) 
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Thus: 






v 

V 

a. 1 



a 



m 



A = 



0 10 0 
0 0 1-1 
0 0 0 0 

ooo -i; 

* 



B = 



O' 

0 

0 

1 

It] 



(Eq 4.68) 




Figure 4.3 Single - lag Guidance System Model 



To find the optimal feedback law from equation 4.65, Y and V have to be calculated using 
equations 4.62 and 4.63, respectively. The solution for Y is obtained from: 

T 



Y + A t Y = 0, 



Y(t F ) = 



0 . 
O’ 
0 



(Eq 4.69) 
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(Eq 4.70) 



0 0 0 0 
10 0 0 

r+oioo^o, nt F ) = 

i 

,0 -1 0 -i\ 

T 



m 

0 

0 

0 



Y is obtained by integrating equation 4.70 backwards, and is: 



Y = 



l F 1 
( t F -t V 



/ x (V-0- 

1 I ( tp — 0 T 

l T J 



go 

'jo 

2 



-r 2 ' e T + 4^ - 1 



(Eq 4.7 1) 



The matrix V (in this case, V is a scalar since we are specifying only one terminal 
constraint, namely zero miss distance: y = 0), is computed using equation 4.63: 

(Eq 4.72) 



V (t) = -jY T BB T Ydt = (T 2 ) 1 1 e r + ^-l 



i 



J 



dt ; where t g0 = t F -t 



After some cumbersome computations, we find: 



V{t) = L (-k’ + Sk 2 -3k) + l2ke~ k + 3e~ 2k -3; 



(Eq 4.73) 

where k = ~ 

T 



Then using equation 4.65 we obtain the optimal feedback control law: 

u(t) = [y + yt go + 0.5a,t 2 + {-a m T 2 ) (e~ k + k- 1)] , (Eq4.74) 



go 



where N = 



6k 2 ( e~ k + k - 1) 



2k 3 -6k 2 + 6k + 3- nke~ k -3e~ 2lc ’ 



k=f- . (Eq 4.75) 



It is desirable to express the state variables y and y in terms of the line of sight rate k . 
Assuming small angles, we fmd from Figure 4.2 that: 
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(Eq 4.76) 





hence, the optimal control law can be expressed as: 



u(r) = NV c X + ^(e k + k- l)a m + ja r 



(Eq 4.77) 



This equation is a biased proportional navigation guidance law where a time varying 
navigation gain N and an acceleration feedback path provide compensation for the missile 
time lag. The acceleration command is issued normal to the LOS. 



In this section we are going to model the missile/target tridimensional engagement 
scenario. Three guidance laws, namely: PROPNAV, augmented PROPNAV and optimal 
guidance will be used and tested in missile guidance. The engagement for the two later' 
guidance laws will be modeled, by assuming the presence of a seeker and a camera aboard 
the missile to extract the target’s 3-D acceleration. The three dimensional MATLAB 
programs are presented in Appendices A through C. The simulation results are presented 
in the next chapter. 

1. 3-D Missile /Target Geometry 

The tridimensional scenario is developed in spherical coordinates by defming two 
perpendicular planes in pitch and yaw, as illustrated by Figure 4.4. 

In Figure 4.4 r is the relative distance between missile and target; ^ pitch and \ aw 

are the line of sight angles over the pitch and yaw planes, respectively, and may be 
expressed as: 



B. TRIDIMENSIONAL MISSILE/TARGET ENGAGEMENT 




( z t~ z m > 



(Eq 4.78) 
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atan 






yaw 



r (y,~_ yJ l 

(v,-.v m ) j' 



(Eq 4.79) 



The coordinate system shown in Figure 4.4 translates with the missile. To track 
the missile target tridimensional positions we define a ground based coordinate system 
shown in Figure 4.5. 




In Figure 4.5: 



X ... = atan 

pitch 



l X m + y m ' 



(Eq 4.80) 
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atan 



(Eq 4.81) 



i_ pitch 



2 2 
■i x J+yjJ 




X 



F igure 4.5 Ground Coordinate System 



X = atan 

m_ yaw 



ov 



(Eq 4.82) 






(Eq 4.83) 



The missile is controlled in 3-D space by issuing guidance commands in two 
orthogonal planes (pitch and yaw), u pi(ch and u yaw . The magnitude of these commands 
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depends on the selected guidance law and their tridimensional direction is perpendicular to 
the LOS defined in the two planes (see Figure 4.4). The LOS in pitch is defined by the 
imaginary 3-D line from the missile to the target in the pitch plane. The LOS in yaw is 
defined by the imaginary 2-D line from the missile to the projection of the target over the 
yaw plane. The yaw plane is simply the horizontal xv plane. The pitch plane is the vertical 
plane normal to the horizontal plane and rotated by the yaw angle X . 

The guidance laws under study can be expressed as a function of the classical 
proportional navigation guidance law plus a term that may depend, among other variables, 
on the target and missile accelerations. Hence, each guidance law (classical proportional 
navigation, augmented proportional navigation and optimal guidance) can, in general, be 
expressed, as: 

u m ( t ) = NV C X +f(a m , a r t g0 , T ) , (Eq 4.84) 

where u m ( t ) is the guidance command that is issued perpendicular to either the LOS in 
yaw or the LOS in pitch. N is constant for PROPNAV and augmented PROPNAV. For 
optimal guidance, N is function of the time to go and the effective time constant of the 
guidance system. The closing speed V c , is the relative speed between the target and the 

missile along either the LOS in yaw or the LOS in pitch. The LOS rate X may be the LOS 
rate in either the pitch or the yaw planes. The term f(a m , a t , t g0 , T) may be function of both 
the missile’s acceleration a m and the target’s acceleration a r the time to go t g0 and the 
guidance system’s effective time constant T. In order to generate the pitch and yaw 
guidance commands u pitch {t) and u yaw (t ) , it is first necessary to explain how to obtain 
the variables that they depend on. 

2. Seeker Head Modeling 

The seeker is able to detect, acquire and track by sensing and processing the 
radiation or reflection of energy by the target. The seeker is normally located in the 
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missile’s nose and mounted on a gimballed platform which maintains the target within the 
field of view by rotating the platform. 

The control torque to the seeker may be described by the following equation. 

r = /p, (Eq 4.85) 

where T is the applied torque, I is the seeker’s moment of inertia and P is the seeker’s 
angular acceleration. The seeker’s dynamics is modeled by the following second order 
differential equation: 

(3 = y = -Cj(P-X) -c 2 P, (Eq 4.86) 

where the coefficients and c 2 are determined by the seeker’s time constant (t .) and 

damping ratio. Taking the Laplace transform of equation 4.86, assuming zero initial 
conditions, we obtain the filter’s transfer function that represents the relationship between 
the LOS angle input A. (s) and the seeker head angle output (3 ( 5 ) : 



P (*) 



^ ( 5 ) s 2 + C 2 S + 



(Eq 4.87) 



Assuming a damping ratio of one, the transfer function in equation 4.87 may be rewritten 
as: 

(3(s) c l c l 



A (s) 



S + C 2 S + C J 



1 2 ' 

^sk 



(Eq 4.88) 



Choosing x sk = 0.1 sec (which is a good approximation of a real world system), the 
constants and c 2 may be obtained: 



1 2 

Cl = ( — ) = 100, 

^ sk 



c 2 = 2 (f-) = 20. 

'"sk 



(Eq 4.89) 
(Eq 4.90) 
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Given that we are interested in the 3-D missile/target engagement, the seeker must provide 
line of sight rate information in both planes. Hence: 



pitch 


20 


(Eq 4.91) 


^ pitch 


s 2 + 20.v + 100 ' 


(3 

^ yaw 


20 


(Eq 4.92) 


X 

yaw 


s 2 + 20s+ 100 



where & puch and (3 yavv are the seeker’s pitch and yaw angles, respectively. Figures 4.6 and 
4.7 depict the pitch and yaw signal flow graphs of the seeker. 




Figure 4.6 Seeker Head Flow Graph (Pitch) 



From these diagrams, the continuous-time state equations of the form: 

k ~ ^sh^sk B sk ll sk' (^9 4.93) 

may be easily obtained. Selecting the state vector to be: 
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Figure 4.7 Seeker Head Flow Graph (Yaw) 



and the seeker head input as: 



u sk = 



pitch 

X 

yaw 



equation 4.93 becomes: 



1 


0 


1 


0 


0 




~ 0 


0 ~ 


X sk = 


-100 


-20 


0 


0 


X sk + 


100 


0 


0 


0 


0 


1 


0 


0 




0 


0 


-100 


-20 




0 


100 



(Eq 4.95) 



(Eq 4.96) 
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The variables f3 pitch and (3 vavv are estimates of the LOS angle rates X pi[ch and 

X aw , and are available from the second and forth states of .v sk , respectively. The estimates 
of the LOS rate in pitch and yaw permits us to determine the missile command inputs u pi(ch 



3. Guidance System 

In this work, the guidance system dynamics are modeled as a single lag as seen in 
equation 4.66. This equation is repeated here. 

a m (*) 1 



U m ( S ) 



1 + 7V 



(Eq 4.97) 



For the 3-D missile/target engagement the guidance system generates missile commands in 
both planes, pitch and yaw. Hence: 

^m. pitch (s') 1 



“ pitch ( s ) 



a m_ yaw 
^ yaw 



1 +7V 
1 

T+7V 



(Eq 4.98) 
(Eq 4.99) 



where a m pitch and a m yaw are the pitch and yaw missile’s accelerations; u pitch and u vaw 

are the pitch and yaw missile’s acceleration commands. Figures 4.8 and 4.9, show the pitch 
and yaw signal flow graphs, for the missile guidance system. We chose the guidance system 
time constant T to be 1.0 second. 
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Figure 4.8 Guidance System Signal Flow Graph (Pitch) 




From these diagrams, the state equation is easily obtained. Defining the guidance system 
state vector as: 



x 



gs 



r 

a m. pilch 
yaw 



and the guidance system input as: 



(Eq 4.100) 
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(Eq 4.101) 



u 



gs 



the state equation becomes: 



11 pitch 
^yaw\ 



X 



gs 



1 0 
-1 



L 0 lj 



gs- 



(Eq 4.102) 



4. Missile And Target Kinematics 

The missile is controlled in three dimensional space by generating acceleration 
commands in two orthogonal planes. These planes are the pitch and yaw planes. The 
acceleration commands in pitch and yaw are issued perpendicular to the respective lines of 
sight. The magnitude of the acceleration commands depends on the selected guidance law. 
From equation 4.84, the pitch and yaw missile acceleration commands may be expressed, 
in its general form, as: 

ll pitch ~NV c _pitch^'pitch'*'fpitch( a m’ a t , tgo , 'r)’ (Eq 4.103) 

u yaw = NV c_yaw\a» + fyaw( a m> a '’ t go> r > ■ (Eq 4.104) 

F c pitch and V c vaw are the relative speeds, between the target and the missile along the 
pitch and yaw line of sights. The functions f pitch ( a m , a (> t g0 , T) and f yaw (a ffl , a ( , t go , T) 

are the augmented PROPNAV or optimal guidance extra terms in pitch and yaw, 
respectively. 

In order to track the missile’s 3-D coordinates {x m , y m , z m ) , the missile command 

accelerations in pitch and yaw, are broken down into cartesian coordinate system 
components. 

Figure 4. 10, shows the decomposition of the pitch acceleration command in its 
components. From this figure the following relationships are derived: 

a mx. pitch ~ ~ ^ a m. pitch^^^pitch^ C0 ^^'yaw' 

^mx. pitch ~ ~ pitch 5 ^^“pitch) 
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(Eq 4.105) 
(Eq 4.106) 



(Eq 4.107) 



a 



mz _ pitch 



a m_ pitch C0S ^- pitch' 




Figure 4.10 Pitch Plane Acceleration Components 



Figure 4.11 shows the decomposition of the yaw plane acceleration command in 
to its x and y components. From the figure the following relationships are obtained: 

(Eq 4.108) 

(Eq 4.109) 



a mx_yaw \vavv’ 



a my_yaw a m_ yaw^^^yaw' 
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Figure 4.11 Yaw Plane Acceleration Components 



To find the overall missile’s acceleration components along the three cartesian axis, we use 
the results of the equations 4.105 through 4.109: 



a mx_ pitch + a mx_ yaw * 


(Eq 4.110) 


^ my _ pitch a my_ yaw ’ 


(Eq 4.111) 


mz ~~ a mz_ pitch’ 


(Eq 4.112) 
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The missile’s tridimensional movement is determined by these three acceleration 
components. Defining the missile state vector as: 

r i 

X m = H (Eq 4.1 13) 

y m 

7 I 

\ m 

L Z ™ 

and the input as the missile’s acceleration command: 






a m - a 



my 



a r 



the missile state equation is: 
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m- 



(Eq 4.114) 



(Eq 4.115) 



As we have seen in Chapter EH, the missile when equipped with a camera and a 
seeker is able to estimate the target’s 3-D acceleration components. This information may 
be used to improve the missile guidance towards the target. Defining the target state vector 
as: 



X 



t 



y t 

y 

y, 

AJ 



(Eq 4.116) 



and the target’s acceleration as: 
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r 



a 



tx 



a t a tv ’ 



the target state equation is: 
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(Eq 4.117) 



(Eq 4.118) 



the tridimensional target acceleration a t may be estimated using the missile’s image 
processing capabilities; 

5. Pitch And Yaw Closing Velocities. Determination Of Time To Go 

Figure 4.12 shows the decomposition of the missile and target absolute velocities. 
From this figure: 



mx. 



V = V 

m r my ’ 

V 



(Eq 4.1 19) 



«J 



[V 

I 1 

v , = IK 



iy 

v 

L '-u 



(Eq 4.120) 



Figures 4. 13 and 4.14 show the projection of the missile’s velocity vector over the pitch and 
yaw planes. 

These figures permit us to compute the missile’s and target’s velocity 
components, over the pitch and yaw planes: 

V m.pitch = ^ cos (y m _ yaw - X yaJ * (Eq 4.121) 



yaw 



= ^COSY 



m ver y 



(Eq 4.122) 



61 



and 



V t_ pitch = V 0 S (Y t.ya W -\aJ> 


(Eq 4.123) 


V t.yaw = V t COS Y t _ ver' 


(Eq 4.124) 




Figure 4.12 Missile And Target Velocity Components 



The pitch closing velocity V c pitch is found by projecting the missile’s and target’s pitch 
plane velocities along the LOS in pitch (see Figure 4.4). Then: 

(Eq 4.125) 

V c. pitch = ^m. pitch C0S ( ^ pitch ~ ^m_ pitch ^ ~ . pitch C0S ( ^ pitch ~ \ pitci) 
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Similarly, the yaw plane closing velocity, is obtained as: 



(Eq 4.126) 



= V 

vow m_ \aw 



C0S ( V yaw ~ KaJ ~ V L yaw* os ( Y,_ vaw ~ KaJ ' 




Figure 4.13 Missile’s Pitch And Yaw Velocity Components 



The time to go until interception may be computed from: 

(Eq 4.127) 

v c_ pitch 

where r is the missile/target relative distance along the pitch LOS (see Figure 4.4). 
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Figure 4.14 Target’s Pitch And Yaw Velocity Components 
6. Proportional Navigation 

The PROPNAV 3-D simulation, in Appendix A, uses missile commands in pitch 
and yaw of the form: 



1 pitch Wc_ pitch^" pitch ’ 


(Eq 4.128) 


l *yaw NV C _ yaw-yaw 


(Eq 4.129) 



where A is a constant. 

7. Augmented Proportional Navigation 

The augmented PROPNAV 3-D simulation, in Appendix B, uses missile 
commands in pitch and yaw of the form: 

11 pitch = NV c. pitch^pitch + 0.5Na { pitch , (Eq 4.130) 
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(Eq 4.131) 



yaw 



li yaw Me. yaw-yaw 0-5 vu 

where a t pltch and a ( yaw are the components of the target ‘s acceleration normal to the 
pitch and yaw line of sights and /V is a constant. 



8. Optimal Guidance 

The optimal guidance 3-D simulation, in Appendix C, uses missile commands in 
pitch and yaw of the form: 



U pitch ^^c.pitch^pitch^ j2^ e ^ a m_ pitch + 2 pitch' (Eq 4.132) 

^yaw ~ c. yaw-yaw "** ^2^ + k ~ 1) yaw + &[_ y aw (Eq 4.133) 

where k and N are given by equation 4.75. 



9. Discrete-Time Simulation Using State Space Methods 
The general continuous-time state equations are: 

x(t) = Ax(t) +Bu(t), (Eq 4.134) 

y(t) = Cx(t ) +Du(t) (Eq 4.135) 

where x (/) is the state vector and y (t) is the output vector. This system is simulated by 
iterating the discrete-time state equations: 

jc(/i+1) = Ojc(«) +ru(n), (Eq 4.136) 

y(n) = Cx(n ) +Du(n). (Eq 4.137) 

where: 



0 > = 




T s = sampling time; 



T, 

A = j 'e A ‘B dt. 
0 



(Eq 4.138) 
(Eq 4.139) 
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The missile/target engagement scenario is simulated 


using the MATLAB 


software package. The discrete-time state equations used in the simulation are defined for 
the seeker, guidance system, missile, and target dynamics: 


x sk (n+l) = * sk x sk (n) + T sk u sk (n); 


(Eq 4.140) 


x gs (n + l) = * gi x gs {n) + r g u gs (n ) ; 


(Eq 4.141) 


.r m (n + l) = d> m .v m (Ai) +T m u m (n); 


(Eq 4.142) 


x t (n+l) = O f .v, (n) + f t u t (/i) . 


(Eq 4.143) 
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V. SIMULATION RESULTS 



A. GENERAL 

This chapter presents the results of the computer simulations for a missile equipped 
with a radar and a camera. The missile/target 3-D engagement was simulated using three 
control laws: 

1 . Proportional navigation, 

2. Augmented proportional navigation, 

3. Optimal guidance. 

The simulation was conducted for two different target maneuvers: 

1. A 3-D constant target acceleration, 

2. A 3-D varying target acceleration. 

The following assumptions are made throughout: 

1. We assume that the simulation’s initial conditions are defined when the missile 

enters into the terminal phase of the flight (about 10 seconds before impact). 

2. The PROPNAV and APROPNAV effective navigation constant is 3 (A = 3), 

3. The missile is limited to 25 g’s accelerations in pitch and yaw, 

4. The instantaneous target acceleration is available from previous image 
processing. No delays are assumed in this process. 

5. The missile and target speeds are limited to 3500 and 2000 feet/sec, 
respectively, 

6. The missile is in a collision triangle with the target on entering into the terminal 
phase of flight. 

7. The target may start its evasion maneuver at any time (initial time), 

8. The target is limited to a maximum of 12 g’s, 

9. The acceleration due to gravity is ignored, 

10. The sampling time is 0.01 seconds. 
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B. ENGAGEMENT SCENARIOS. RESULTS 



1. Scenario #1 (Constant Target Acceleration) 

The initial missile/target geometry (when the missile enters into the terminal 
phase of flight) is shown in Figure 5.1. 




Figure 5.1 Initial Missile/Target Geometry 



The engagement initial conditions are (distances are in feet and speeds in feet / 

sec): 
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(Eq 5.1) 



(Eq 5.2) 



The target’s evasive maneuver is constant (the accelerations are in feet/sec A 2) 





3 xgj 


y, 


= 1 xg 


-L 


L 2 *& 



(Eq 5.3) 



'j 

where the acceleration of gravity is: g = 32.2 feet/sec . Figures 5.2 to 5.22 display the 
results of the three dimensional simulation for the constant target evasive maneuver. 
Figures 5.2 to 5.8 relate to the proportional navigation (PROPNAV) guidance law. Figures 
5.9 to 5.15 relate to the augmented proportional (APROPNAV) guidance law. Figures 5.16 
to 5.22 relate to the optimal guidance law. Figures 5.4 to 5.8, 5.1 1 to 5.15 and 5.18 to 5.22 
display the results assuming that the target starts its maneuver 6 seconds after the missile 
entered into the terminal phase of flight. 



69 





70 



MISSILE ACCELERATION MAGNITUDE vs TIME, PROPNAV 




Figure 5.4 Missile Acceleration Magnitude (PROPNAV) 




Figure 5.5 Target Acceleration Magnitude (PROPNAV) 
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Figure 5.7 Missile Yaw Acceleration (PROPNAV) 
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Figure 5.9 Miss Distance vs. Time To Go (APROPNAV) 



MISS DISTANCE vs INITIAL TIME, APROPNAV 




Figure 5.10 Miss Distance vs. Initial Time (APROPNAV) 
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Figure 5.12 Target Acceleration Magnitude (APROPNAV) 
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Figure 5.13 Missile Pitch Acceleration (APROPNAV) 




Figure 5.14 Missile Yaw Acceleration (APROPNAV) 
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MISS DISTANCE vs TIME TO GO. OPTIMAL 




Figure 5.16 Miss Distance vs. Time To Go (OPTIMAL) 



MISS DISTANCE vs INITIAL TIME. OPTIMAL 




Figure 5.17 Miss Distance vs. Initial Time (OPTIMAL) 
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MISSILE PITCH ACCELERATION vs TIME. OPTIMAL 




Figure 5.20 Missile Pitch Acceleration (OPTIMAL) 




Figure 5.21 Missile Yaw Acceleration (OPTIMAL) 
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2. Scenario #2 (Varying Target Acceleration) 

The initial missile/target geometry and the initial conditions are the same as used 
for the constant acceleration scenario. Namely, the target geometry is shown in Figure 5.1 
and the initial conditions are defined by equations 5.4 and 5.5. The target’s evasive 
maneuver may start at any time to go. The target performs the following 3-D maneuver (the 
accelerations are in feet/sec A 2): 



where y ( vaw and y t pitch are the target’s yaw and pitch flight path angles, respectively. 

Figures 5.23 to 5.43 display the results of the three dimensional simulation for this type of 
evasive maneuver. Figures 5.23 to 5.29 relate to the proportional navigation (PROPNAV) 
guidance law. Figures 5.30 to 5.36 relate to the augmented proportional (APROPNAV) 
guidance law. Figures 5.37 to 5.43 relate to the optimal guidance law. Figures 5.25 to 5.29, 
5.32 to 5.36 and 5.39 to 5.43 display results, assuming that the target starts its maneuver 6 
seconds after the missile entered into the terminal phase of flight. 



(Eq 5.4) 




.V, = 4 x < gxcos(Y f _ yaw ) 
h 3 x g x cos (y t pitch ) 
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Figure 5.25 Missile Acceleration Magnitude (PROPNAV) 




Figure 5.26 Target Acceleration Magnitude (PROPNAV) 
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Figure 5.29 Missile To Target Range (PROPNAV) 
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Figure 5.31 Miss Distance vs. Initial Time (APROPNAV) 
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MISSILE ACCELERATION MAGNITUDE vs TIME. APROPNAV 




Figure 5.32 Missile Acceleration Magnitude (APROPNAV) 




Figure 5.33 Target Acceleration Magnitude (APROPNAV) 
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Figure 5.34 Missile Pitch Acceleration (APROPNAV) 




Figure 5.35 Missile Yaw Acceleration (APROPNAV) 
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MISS DISTANCE vs TIME TO GO. OPTIMAL 




Figure 5.37 Miss Distance vs. Time To Go (OPTIMAL) 



MISS DISTANCE vs INITIAL TIME, OPTIMAL 




Figure5.38 Miss Distance vs. Initial Time (OPTIMAL) 
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Figure 5.39 Missile Acceleration Magnitude (OPTIMAL) 
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TARGET ACCELERATION MAGNITUDE vs TIME. OPTIMAL 
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Figure 5.40 Target Acceleration Magnitude (OPTIMAL) 
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Figure 5.41 Missile Pitch Acceleration (OPTIMAL) 




Figure 5.42 Missile Yaw Acceleration (OPTIMAL) 
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VI. CONCLUSIONS AND RECOMENDATIONS 



A. CONCLUSIONS 

l. Scenario #1 (Constant Target Acceleration) 




Figure 5.44 Miss Distance Comparison For The Three Guidance Laws 

(Scenario #1) 
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Figure 5.45 Missile Acceleration Comparison For The Three Guidance Laws 

(Scenario #1) 



% 



2. Scenario #2(Varying Target Acceleration) 




Figure 5.46 Miss Distance Comparison For The Three Guidance Laws 

(Scenario #2) 
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Figure 5.47 Missile Acceleration Comparison For The Three Guidance Laws 

(Scenario #2) 
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3. Discussion 

The three dimensional miss distance may be improved by estimating the 3-D 
target’s evasive maneuver. One way to estimate the 3-D target acceleration is by utilizing 
dynamic image processing. Three setups were considered: 

1 . two cameras, 

2. a camera and a radar, 

3. only one camera. 

In the single camera setup, actual range is achieved by guessing the target’s 
physical dimensions. The target’s 3-D motion parameters can be estimated by utilizing two 
consecutive image frames. The target’s acceleration may be computed by taking the second 
derivative after filtering the target’s motion data. Alternatively, the target’s 3-D motion 
may be processed by Kalman filters to estimate its acceleration components. This 
additional information about the target’s behavior is injected in suitable control laws to 
improve the missile’s homing performance. 

Proportional navigation, augmented proportional navigation, and optimal 
guidance laws were derived for use in a three dimensional environment. The classical 
proportional navigation guidance law tracks a target with good accuracy, especially if the 
target maneuvers at long time to go. However, when compared with augmented 
PROPNAV and optimal guidance, PROPNAV requires higher missile acceleration 
capabilities. A plausible guidance law is one that issues missile’s commands proportional 
to the miss distance that would result if the missile made no further corrections. Augmented 
proportional navigation was derived using this heuristic argument. For a constant target 
maneuver, augmented proportional navigation increases the missile percentage of kill. For 
a non constant evasive maneuver, APROPNAV does not always guarantees less miss 
distance than PROPNAV. However, APROPNAV requires less missile acceleration 
capabilities than PROPNAV. Optimal guidance was derived for a missile with a single lag 
guidance system. Optimal guidance provides compensation for the missile’s guidance 
system dynamics. The optimal guidance law requires the least missile acceleration 
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capability of the three guidance laws. In fact, this law is derived in order to drive the miss 
distance to zero while minimizing a performance index made up of the integral of the 
square of the control. Clearly, the optimal guidance law presents the least miss distance of 
the three guidance laws. However, it requires a missile with complex signal processing 
capabilities. The homing capabilities of the missile can be dramatically increased by 
identifying the target’s evasive maneuver and injecting this information into the 
APROPNAV (especially for a constant target maneuver) or optimal guidance control 
algorithms. The optimal control algorithm guarantees extraordinary performance. Utilizing 
optimal guidance, especially against highly responsive targets, can be the difference 
between failure and success. 



B. RECOMEND ATIONS 

It is recommended to continue this research by simulating the overall system (i.e. 
estimating the 3-D target’s evasive maneuver from two consecutive image frames and 
injecting this data into the tridimensional missile/target engagement simulation programs 
developed in this thesis). The simulations developed in this thesis are very generic and 
easily adapted to different conditions (i.e. for systems with different dynamics and initial 
conditions). The consequences of the image measurement errors in the target acceleration 
estimation and ultimately in the miss distance can be investigated. Finally, it is 
recommended that electronic counter measures (ECCM) be added to the target’s evasion 
maneuver in order to evaluate their effects on the miss distance. 
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APPENDIX A - MISSILE/TARGET THREE DIMENSIONAL 
SIMULATION USING PROPORTIONAL NAVIGATION GUIDANCE 



% Written by: Rui Manuel Alves Francisco 
% Date: 14 October 1992 

% This Program simulates the terminal phase of a 3-D missile/target 
% engagement using classical proportional navigation. 

clear 

clg 

% DEFINE CONSTANTS 
dt = .01; % Sampling time 
Tf = 100; % maximum simulation time 
kmax = Tf/dt+1; 
n = 3; % navigation constant 
N = [nO 
On]; 

% DEFINE STATE EQUATIONS 
% Missile 

% Xm = [xm = missile ‘s x coordinate 
% xmd = missile's speed (x coordinate) 

% ym = missile's y coordinate 

% ymd = missile's speed (y coordinate) 

% zm = missile's z coordinate 

% zmd = missile's speed (z coordinate)] 

Am = [0 1 0 0 0 0 
000000 
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000100 

000000 

000001 

000000 ]; 

Bm = [0 0 0 
100 
000 
0 1 0 
000 
0 0 1 ]; 

% Target 

% Xt = [xt = target’s x coordinate 
% xtd = target’s speed (x coordinate) 
% yt = target’s y coordinate 
% ytd = target’s speed (y coordinate) 
% zt = target’s z coordinate 
% ztd = target’s speed (z coordinate)] 

At = [0 1 0 0 0 0 
000000 
000100 
000000 
000001 
00 0000 ]; 

Bt = [0 0 0 
100 
000 
0 1 0 
000 
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0 0 1 ]; 

% Seeker (Radar) 

% Xsk = [beta_pitch = seeker’s pitch angle 
% betad_pitch = seeker’s pitch angle rate 

% beta_yaw = seeker’s yaw angle 

% betad_yaw = seeker’s yaw angle rate] 

Ask = [ 0 1 0 0 

-100 -20 0 0 

0 0 0 1 

0 0 -100 -20]; 

Bsk = [00 
100 0 
0 0 
0 100 ]; 

% Guidance System 

% Xgs = [a_m_pitch = missile’s pitch acceleration 
% a_m_yaw = missile’s yaw acceleration] 

Ags = [-10 

0-1]; 

Bgs = [10 

0 1 ]; 

% INITIALIZE STATE VARIABLES (when the missile enters into the terminal 
phase of flight) 

% Missile 

Xm(:,l) = ( 0 % The missile is in a collision triangle 

2828 % with the target when the missile enters into 

0 % the terminal phase of flight 

1000 
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0 

47.1339]; 

% Target 
Xt(:,l) = [30000 
0 
0 

1000 

500 

0 ]; 

% DISCRETE REPRESENTATION 
[PHIm.DELm] = c2d(Am,Bm,dt); 

[PHIt.DELt] = c2d(At,Bt,dt); 

[PHIsk.DELsk] = c2d(Ask,Bsk,dt); 

[PHIgs.DELgs] = c2d(Ags,Bgs,dt); 

%LINE OF SIGHT (LOS) INFORMATION. INITIAL CONDITIONS. 
% Missile 

% LAMBDA_m = Missile’s LOS from the global coordinate system 
% LAMBDA_m = [LAMB DA_m_p itch = Missile’s pitch LOS angle 
% LAMB D A_m_y aw = Missile’s yaw LOS angle] 

LAMBDA_m(:,l) = [atan2(Xm(5,l),sqrt(Xm(l,l) A 2+Xm(3,l) A 2)); 

atan2(Xm(3, 1 ),Xm( 1,1))]; 

% Target 

% LAMBDA_t = Target’s LOS from the global coordinate system 
% LAMBDA_t = [LAMBDA_t_pitch = Target’s pitch LOS angle 
% LAMBDA_t_yaw = Target’s yaw LOS angle] 

LAMBDA_t(:,l) = [atan2(Xt(5,l),sqrt(Xt(l,l) A 2+Xt(3,l) A 2)); 

atan2(Xt(3,l),Xt(l,l))]; 

% LOS from Missile to Target 
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% LAMBDA = LOS from Missile to Target. 

% LAMBDA = [LAMBDA_pitch = LOS angle in pitch 
% LAMBDA_yaw = LOS angle in yaw]; 

LAMBDA(:,1) = [atan2((Xt(5,l)-Xm(5,l)),sqrt((Xt(l,l)-Xm(l,l)) A 2 
+(Xt(3,l)-Xm(3,l)) A 2)); 

atan2((Xt(3, 1 )-Xm(3, 1 )) ,abs(Xt( 1 , 1 )-Xm( 1 , 1 )))] ; 

% MISSILE and TARGET FLIGHT PATH ANGLES INFORMATION 
% Missile 

%GAMMA_m = [GAMMA_m_pitch = Missile’s flight path angle in pitch 
% GAMMA_m_yaw = Missile’s flight path angle in yaw] 

GAMMA_m(:,l) = [atan2(Xm(6,l),sqrt(Xm(2,l) A 2+Xm(4,l) A 2)); 
atan2(Xm(4,l),Xm(2,l))]; 

% Target 

%GAMMA_t = [GAMMA_t_pitch = Target’s flight path angle in pitch 
% GAMMA_t_YAW = Target’s flight path angle in yaw] 

GAMMA_t(:,l) = [atan2(Xt(6,l),sqrt(Xt(2,l) A 2+Xt(4,l) A 2)); 

atan2 (Xt(4, 1 ) ,Xt(2 , 1 ))] ; 

% RANGE INFORMATION 
% Missile 

% Rm = Missile’s range 

Rm(l) = sqrt(Xm(l,l) A 2 + Xm(3,l) A 2 + Xm(5,l) A 2); 

% Target 

% Rt = Target’s range 

Rt(l) = sqrt(Xt(l,l) A 2 + Xt(3,l) A 2 + Xt(5,l) A 2); 

% Missile/Target relative distance 
% R = [Rmtx = Missile/Target x coordinate range 
% Rmty = Missile/Target y coordinate range 

% Rmtz = Missile/Target z coordinate range 
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% Rmt = Missile/Target relative distance(Rmt=sqrt(Rmtx A 2+Rmty A 2+Rmtz A 2))] 
R(:,l) = [Xt(l,l)-Xm(l,l) 

Xt(3,l)-Xm(3,l) 

Xt(5,l)-Xm(5,l) 

sqrt((Xt(l,l)-Xm(l,l)) A 2+(Xt(3,l)-Xm(3,l)) A 2+(Xt(5,l)-Xm(5,l)) A 2)]; 

% VELOCITY INFORMATION 
% Missile 

% Vm = Missile’s Speed 

Vm(l) = sqrt(Xm(2,I) A 2+Xm(4,l) A 2+Xm(6,l) A 2); 

% Target 

% Vt = Target’s Speed 

Vt(l) = sqrt(Xt(2, 1 ) A 2+Xt(4, 1 ) A 2+Xt(6, 1 ) A 2); 

% Speed along the pitch and yaw LOS. Pitch and yaw closing speeds 
Vt_pitch(l) = Vt ( 1 ) *cos (LAMB D A( 2 , 1 )-G AMMA_t(2 , 1 ) ) ; 

Vm_pitch(l) = Vm(l)*cos(LAMBDA(2,l)-GAMMA_m(2,l)); 

Vc_pitch(l) = Vm_pitch( 1 )*cos(GAMMA_m( 1 , 1)-LAMBDA( 1,1)) 
-Vt_pitch(l)*cos(GAMMA_t(l,l)-LAMBDA(l,l)); 

Vt_yaw(l) = V t( 1 ) *cos(G AMMA_t( 1,1)); 

Vm_yaw(l) = Vm(l)*cos(GAMMA_m(l,l)); 

Vc_yaw(l) = Vm_yaw(l)*cos(GAMMA_m(2,l)-LAMBDA(2,l)) 

- Vt_yaw( 1 )*cos(GAMMA_t(2, 1 )-LAMBDA(2, 1 )); 

Vc = [Vc_pitch(l) 0 

0 Vc_yaw(l)]; 

% SEEKER and GUIDANCE SYSTEM INITIAL CONDITIONS and INPUTS. 
% Seeker 
Xsk(:,l) = [0 
0 
0 



106 



0 ]; 



Usk(:,l) = LAMBDA(:,1); % Seeker input 
% Guidance System 
Xgs(:,l) = [0 
0 ]; 

Ugs(:,l) = N*Vc*[Xsk(2,l) % Guidance system input 
Xsk(4, 1)]; 

% TIME 

% Time = Time vector 
TIME(l) = 0; 

% Tgo = Time to go 
Tgo(l) = R(4,l)/Vc_pitch(l); 

% SIMULATE THE SYSTEM 
for ti = 0:. 25:21.25 
for i = I:kmax-1 

% Calculate components of the missile’s pitch acceleration 
a_mx_pitch(i) = -(Xgs(l,i)*sin(LAMBDA(l,i))*cos(LAMBDA(2,i))); 
a_my_pitch(i) = -(Xgs(l,i)*sin(LAMBDA(l,i))*sin(LAMBDA(2,i))); 
a_mz_pitch(i) = Xgs(l,i)*cos(LAMBDA(l,i)); 

% Calculate components of the missile’s yaw acceleration 
a_mx_yaw(i) = -Xgs(2,i)*sin(LAMBDA(2,i)); 
a_my_yaw(i) = Xgs(2,i)*cos(LAMBDA(2,i)); 

% Compute overall missile acceleration 
a_mx(i) = a_mx_pitch(i)+a_mx_yaw(i); 
a_my(i) = a_my_pitch(i)+a_my_yaw(i); 
a_mz(i) = a_mz_pitch(i); 
am(:,i) = [a_mx(i) 
a_my(i) 
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a_mz(i)]; 

% Compute missile’s acceleration magnitude 
a_m(i) = sqrt(a_mx(i) A 2 + a_my(i) A 2 + a_mz(i) A 2); 

% Generate target’s evasive maneuver (we assume that these accelerations, along 
% the three cartesian axis, are estimated using the missile’s image processing 
% capabilities) 

if TIME(i) >= ti/2 % target starts evasive maneuver 
a_tx(i) = 3*32.2*sin(GAMMA_t(2,i)); 
a_ty(i) = 4*32.2*cos(GAMMA_t(2,i)); 
a_tz(i) = 3*32.2*cos(GAMMA_t(l,i)); 
else 

a_tx(i) = 0.0; 
a_ty(i) = 0.0; 
a_tz(i) = 0.0; 
end 

at(:,i) = [a_tx(i) 
a_ty(i) 
a_tz(i)]; 

% Compute magnitude of the target’s acceleration 
a_t(i) = sqrt(a_tx(i) A 2 + a_ty(i) A 2 + a_tz(i) A 2); 

% Update missile states 

Xm(:,i+1) = PHIm*Xm(;,i)+DELm*am(:,i); 

% Update target states 

Xt(:,i+1) = PHIt*Xt(:,i)+DELt*at(:,i); 

% Update seeker states 

Xsk(:,i+1) = PHIsk*Xsk(:,i)+DELsk*Usk(:,i); 

% Update Guidance System states 
Xgs(:,i+1) = PHIgs*Xgs(:,i)+DELgs*Ugs(:,i); 
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% Limit yaw and pitch accelerations to 25 g’s 
if abs(Xgs(l,i+l)) > 805.0 
Xgs(l,i+1) = 805.0 *sign(Xgs(l,i+l)); 

end 

if abs(Xgs(2,i+l)) >805.0 
Xgs(2,i+1) = 805.0 *sign(Xgs(2,i+l)); 
end 

% Update LOS angles 

LAMBDA(:,i+l) = [atan2((Xt(5,i+l)-Xm(5,i+l)),sqrt((Xm(l,i+l)-Xt(l,i+l)) A 2 
+(Xm(3,i+l)-Xt(3,i+l)) A 2)); 

atan2((Xt(3,i+l)-Xm(3,i+l)),(abs(Xm(l,i+l)-Xt(l,i+l))))]; 
Usk(:,i+1) = LAMBDA(:,i+l); 

LAMBDA_m(:,i+l) = [atan2(Xm(5,i+l),sqrt(Xm(l,i+l) A 2+Xm(3,i+l) A 2)); 

atan2(Xm(3,i+ 1 ),Xm( 1 ,i+ 1))] ; 

LAMBDA_t(:,i+l) = [atan2(Xt(5,i+l),sqrt(Xt(l,i+l) A 2+Xt(3,i+l) A 2)); 

atan2(Xt(3,i+l),Xt(l,i+l))]; 

% Update flight path angles 

GAMMA_m(:,i+l) = [atan2(Xm(6,i+l),sqrt(Xm(2,i+l) A 2+Xm(4,i+l) A 2)); 

atan2(Xm(4,i+ l),Xm(2,i+ 1 ))] ; 

GAMMA_t(.\i+l) = [atan2(Xt(6,i+l),sqrt(Xt(2,i+l) A 2+Xt(4,i+l) A 2)); 

atan2(Xt(4,i+ 1 ) ,Xt(2,i+ 1 ))] ; 

% Update Range Information 

Rm(i+1) = sqrt(Xm(l,i+l) A 2 + Xm(3,i+1) A 2 + Xm(5,i+1) A 2); 

Rt(i+1) = sqrt(Xt(l,i+l) A 2 + Xt(3,i+1) A 2 + Xt(5,i+1) A 2); 

R(:,i+1) = [Xt(l,i+l)-Xm(U+l); 

Xt(3,i+1)-Xm(3,i+1); 

Xt(5 ,i+ 1 ) -Xm (5 ,i+ 1 ) ; 

sqr t((Xt( 14+ 1)-Xm( 1 ,i+ 1 )) A 2+(Xt(3 ,i+ 1 )-Xm(3 ,i+ 1 )) A 2 
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+(Xt(5,i+l)-Xm(5,i+l)) A 2)]; 

% Update Velocity Information 

Vm(i+1) = sqrt(Xm(2,i+ 1 ) A 2+Xm(4,i+ 1 ) A 2+Xm(6,i+ 1 ) A 2); 

Vt(i+1) = sqrt(Xt(2,i+ 1 ) A 2+Xt(4,i+ 1 ) A 2+Xt(6,i+ 1) A 2); 

Vt_pitch(i+1) = Vt(i+l)*cos(LAMBDA(2,i+l)-GAMMA_t(2,i+l)); 
Vm_pitch(i+1) = Vm(i+l)*cos(LAMBDA(2,i+l)-GAMMA_m(2,i+l)); 
Vc_pitch(i+ 1 ) = Vm_pitch(i+l)*cos(GAMMA_m(l,i+l)-LAMBDA(l,i+l)) 
- Vt_pitch(i+ 1 ) *cos (G AMMA_t( 1 ,i+ 1 )-LAMB D A( 1 ,i+ 1 )) ; 
Vt_yaw(i+1) = Vt(i+l)*cos(GAMMA_t(l,i+l)); 

Vm_yaw(i+1) = Vm(i+l)*cos(GAMMA_m(l,i+l)); 

Vc_yaw(i+1) = Vm_yaw(i+l)*cos(GAMMA_m(2,i+l)-LAMBDA(2,i+l)) 
-Vt_yaw(i+l)*cos(GAMMA_t(2,i+l)-LAMBDA(2,i+l)); 
Vc = [Vc_pitch(i+1) 0 

0 Vc_yaw(i+1)]; 

% Update guidance system input 
Ugs(:,i+1) = N*Vc*[Xsk(2,i+l) 

Xsk(4,i+1)]; 

% Update Time/time to go 
TEME(i+l) = TTME(i)+dt; 

Tgo(i+l) = R(4,i+ 1 )/V c_pitch(i+ 1 ); 

% Check for closest point 
if (R(4,i) < R(4,i+1)), break, end 
end; 

% Save information for plotting and evaluation 

Rl(4*ti+1) = R(4,i); % miss distance 

Ti(4*ti+1) = ti/2;% starting time of evasive maneuver (EM) 

tgo(4*ti-i-l) = i*dt-ti/2; % time to go until end of flight 

if ti == 12.0 % Record information for a target that 
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% initialized the evasive maneuver 6 sec after the missile 
% entered into the terminal phase of flight. 

TGO = tgo(49); 

Xseeker = Xsk(:,l:i); 

Xgsys = Xgs(:,l:i); 
lambda_m = LAMBDA_m(:,l:i); 
lambda_t = LAMBDA_t(:,l:i); 
lambda = LAMB DA(:,l:i); 
gamma_m = GAMMA_m(;,l:i); 
gamma_t = GAMMA_t(:,l:i); 
r = R(:,l:i); 
vm = Vm(l.i); 
vt = Vt(l:i); 

vm_pitch = Vm_pitch(l:i); 
vt_pitch = Vt_pitch(l:i); 
vm_yaw = Vm_yaw(l:i); 
vt_yaw = Vt_yaw(l:i); 
vc_pitch = Vc_pitch(l:i); 
vc_yaw = Vc_yaw(l:i); 
tGO = Tgo(l:i); 
a_M = am(:,l:i); 
a_T = at(:,l:i); 

A_t = a_t(l:i); 

A_m = a_m(l:i); 
time = TIME(l:i); 
end 

clear R; 

R(:,l) = [Xt(l,l)-Xm(l,l); 
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Xt(3,l)-Xm(3,l); 

Xt(5,l)-Xm(5,l); 

sqrt((Xt(l,l)-Xm(l,l)) A 2+(Xt(3,l)-Xm(3,l)) A 2+(Xt(5,l)-Xm(5,l)) A 2)]; 

end; 

save thesis lp343 R1 tgo Ti tGO missile target TGO Xseeker Xgsys lambda_m 
lambda_t lambda gamma_m gamma_t r vm vt vm_pitch vm_yaw vt_pitch 
vt_yaw vc_pitch vc_yaw a_M a_T time A_t A_m 
PLOTS 

% Miss distance information 

plot(Ti,R 1 ),title( ‘MISS DISTANCE vs INITIAL TIME, PROPNAV’) 

xlabeK'INITIAL TIME - SEC’),ylabel(‘MISS - FEET’) 

print -dps Rlapl 

Ipstoepsi Rlapl.ps Rlapl. epsi 

pause,clg 

plot(tgo,Rl),titIe(‘MISS DISTANCE vs TIME TO GO, PROPNAV’) 

xlabel(‘TIME TO GO - SEC’),ylabel(‘MISS - FEET’) 

print -dps Rlbpl 

Ipstoepsi Rlbpl.ps Rlbpl. epsi 

pause,clg 

% Missile acceleration information 

plot(time,A_m), title( ‘MISSILE ACCELERATION MAGNITUDE vs TIME, 

PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps A_mp 1 
Ipstoepsi A_mpl.ps A_mpl.epsi 
pause, clg 

plot(time,Xgsys(l,:)),title(‘MISSILE PITCH ACCELERATION vs TIME, 

PROPNAV’) 
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xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps Xgsyslpl 
Ipstoepsi Xgsyslpl.ps Xgsyslpl.epsi 
pause, clg 

plot(time,Xgsys(2,:)),title('MISSILE YAW ACCELERATION vs TIME, 

FROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps Xgsys2pl 
Ipstoepsi Xgsys2pl.ps Xgsys2pl.epsi 
pause, clg 

% Target acceleration information 

plot(time,A_t),title( TARGET ACCELERATION MAGNITUDE vs TIME, 
PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps A_tp 1 
I pstoepsi A_tp 1 .ps A_tp 1 .epsi 
pause, clg 

% Seeker pitch and yaw angles 

plot(time,Xseeker( 1 ,:)),title(‘SEEKER PITCH ANGLE vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD’) 

print -dps Xseekerlpl 

Ipstoepsi Xseekerlpl.ps Xseekerlpl. epsi 

pause, clg 

plot(time,Xseeker(3,:)),title(‘SEEKER YAW ANGLE vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD’) 

print -dps Xseeker2p 1 

Ipstoepsi Xseeker2pl.ps Xseeker2pl.epsi 

pause, clg 
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plot(time,Xseeker(2,:)),title(‘SEEKER PITCH ANGLE RATE vs TIME, 

PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD/SEC’) 
print -dps Xseeker3pl 
Ipstoepsi Xseeker3pl.ps Xseeker3pl.epsi 
pause,clg 

plot(time,Xseeker(4,:)),title(‘SEEKER YAW ANGLE RATE vs TIME, 

PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD/SEC’) 
print -dps Xseeker4pl 
Ipstoepsi Xseeker4pl.ps Xseeker4pl.epsi 
pause.clg 

% Range information 

plot(time,r(4,:)),title(‘RANGE vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET’) 

print -dps rpl 

Ipstoepsi rpl.ps rpl.epsi 

pause.clg 

% Missile velocity information 

plot(time,vm),title(‘MISSILE SPEED vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vmpl 

Ipstoepsi vmpl.ps vmpl.epsi 

pause.clg 

% Target velocity information 

plot( time.vt) ,title( ‘T ARGET SPEED vs TIME, PROPNAV’) 
xlabel(‘TTME - SEC’),ylabel(‘FEET/SEC’) 
print -dps vtp 1 
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Ipstoepsi vtpl.ps vtpl.epsi 
pause,clg 

% Closing velocity information 

plot(time,vc_pitch), title( ‘ PITCH CLOSING SPEED, PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vclpl 

Ipstoepsi vclpl.ps vclpl. epsi 

pause.clg 

plot(time,vc_yaw),title(‘YAW CLOSING SPEED vs TIME, PROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vc2pl 

Ipstoepsi vc2pl.ps vc2pl.epsi 

pause.clg 
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APPENDIX B - MISSILE/TARGET THREE DIMENSIONAL 
SIMULATION USING AUGMENTED PROPORTIONAL NAVIGATION 
GUIDANCE 



% Written by: Rui Manuel Alves Francisco 
% Date: 10 November 1992 

% This Program simulates the terminal phase of a 3-D missile/target 
% engagement using augmented proportional navigation guidance. 

clear 

clg 

% DEFINE CONSTANTS 
dt = .01; % Sampling time 
Tf = 100; % maximum simulation time 
kmax = Tf/dt+1; 
n = 3; % navigation constant 
N = [n 0 
On]; 

% DEFINE STATE EQUATIONS 
% Missile 

% Xm = [xm = missile‘s x coordinate 
% xmd = missile's speed (x coordinate) 

% ym = missile's y coordinate 

% ymd = missile's speed (y coordinate) 

% zm = missile's z coordinate 

% zmd = missile's speed (z coordinate)] 

Am = [0 1 0 0 0 0 



116 



000000 
000100 
000000 
000001 
0 0 00 00 ]; 

Bm = [000 
100 
000 
0 1 0 
000 
0 0 1 ]; 

% Target 

% Xt = [xt = target’s x coordinate 
% xtd = target’s speed (x coordinate) 
% yt = target’s y coordinate 
% ytd = target’s speed (y coordinate) 

% zt = target’s z coordinate 

% ztd = target’s speed (z coordinate)] 

At = [0 1 0 0 0 0 
000000 
000 1 00 
000000 
000001 
0 0 0 0 0 0 ]; 

Bt = [0 0 0 
100 
000 
0 1 0 
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000 
0 0 1 ]; 

% Seeker (Radar) 

% Xsk = [beta_pitch = seeker’s pitch angle 
% betad_pitch = seeker’s pitch angle rate 

% beta_yaw = seeker’s yaw angle 

% betad_yaw = seeker’s yaw angle rate] 

Ask = [ 0 1 0 0 

-100 -20 0 0 

0 0 0 1 

0 0 -100 -20]; 

Bsk = [00 
100 0 
0 0 
0 100 ]; 

% Guidance System 

% Xgs = [a_m_pitch = missile’s pitch acceleration 
% a_m_yaw = missile’s yaw acceleration] 

Ags = [-1 0 

0-1]; 

Bgs = [10 
0 1]; 

% INITIALIZE STATE VARIABLES (when the missile enters into the terminal 
phase of flight) 

% Missile 

Xm(:,l) = [ 0 % The missile is in a collision triangle 

2828 % with the target when the missile enters into 

0 % the terminal phase of flight 
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1000 

0 

47.1339]; 

% Target 
Xt(:,l) = [30000 
0 
0 

1000 

500 

0 ]; 

% DISCRETE REPRESENTATION 
[PHIm,DELm] = c2d(Am,Bm,dt); 

[PHIt,DELt] = c2d(At,Bt,dt); 

[PHIsk,DELsk] = c2d(Ask,Bsk,dt); 

[PHIgs,DELgs] = c2d(Ags,Bgs,dt); 

% LINE OF SIGHT (LOS) INFORMATION. INITIAL CONDITIONS. 
% Missile 

% LAMBDA_m = Missile’s LOS from the global coordinate system 
% LAMBDA_m = [LAMB DA_m_p itch = Missile’s pitch LOS angle 
% LAMB D A_m_y aw = Missile’s yaw LOS angle] 

LAMBDA_m(:,l) = [atan2(Xm(5,l),sqrt(Xm(l,l) A 2+Xm(3,l) A 2)); 

atan2(Xm(3,l),Xm(l,l))]; 

% Target 

% LAMBDA_t = Target’s LOS from the global coordinate system 
% LAMBDA_t = [LAMBDA_t_pitch = Target’s pitch LOS angle 
% LAMB DA_t_yaw = Target’s yaw LOS angle] 

LAMBDA_t(:,l) = [atan2(Xt(5,l),sqrt(Xt(U) A 2+Xt(3,l) A 2)); 
atan2(Xt(3 , 1 ) ,Xt( 1,1))]; 
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% LOS from Missile to Target 
% LAMBDA = LOS from Missile to Target. 

% LAMBDA = [LAMB DA_p itch = LOS angle in pitch 
% LAMBDA_yaw = LOS angle in yaw]; 

LAMBDA(:,1) = [atan2((Xt(5,l)-Xm(5,l)),sqrt((Xt(U)-Xm(l,l)) A 2 
+(Xt(3,l)-Xm(3,l)) A 2)); 

atan2((Xt(3J)-Xm(3J)),abs(Xt(U)-Xm(l,l)))]; 

% MISSILE and TARGET FLIGHT PATH ANGLES INFORMATION 
% Missile 

%GAMMA_m = [GAMMA_m_pitch = Missile’s flight path angle in pitch 
% GAMMA_m_yaw = Missile’s flight path angle in yaw] 

GAMMA_m(:,l) = [atan2(Xm(6,l),sqrt(Xm(2,l) A 2+Xm(4,l) A 2)); 
atan2(Xm(4, 1 ),Xm(2, 1 ))] ; 

% Target 

%GAMMA_t = [GAMMA_t_pitch = Target’s flight path angle in pitch 
% G AMMA_t_Y AW = Target’s flight path angle in yaw] 

GAMMA_t(:,l) = [atan2(Xt(6,l),sqrt(Xt(2,l) A 2+Xt(4,l) A 2)); 

atan2(Xt(4,l),Xt(2,l))]; 

% RANGE INFORMATION 
% Missile 

% Rm = Missile’s range 

Rm(l) = sqrt(Xm(l,l) A 2 + Xm(3,l) A 2 + Xm(5,l) A 2); 

% Target 

% Rt = Target’s range 

Rt(l) = sqrt(Xt(l,l) A 2 + Xt(3,l) A 2 + Xt(5,l) A 2); 

% Missile/Target relative distance 
% R = [Rmtx = Missile/Target x coordinate range 
% Rmty = Missile/Target y coordinate range 
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% Rmtz = Missile/Target z coordinate range 

% Rmt = Missile/Target relative distance(Rmt=sqrt(Rmtx A 2+Rmty A 2+Rmtz A 2))] 
R(:,l) = [Xt(l,l)-Xm(l,l) 

Xt(3,l)-Xm(3,l) 

Xt(5,l)-Xm(5,l) 

sqrt((Xt( 1 , 1 )-Xm( 1 , l)) A 2+(Xt(3, 1 )-Xm(3, 1 )) A 2+(Xt(5, 1 )-Xm(5, 1 )) A 2)] ; 

% VELOCITY INFORMATION 
% Missile 

% Vm = Missile’s Speed 

Vm(l) = sqrt(Xm(2 , 1 ) A 2+Xm(4, 1 ) A 2+Xm(6, 1 ) A 2); 

% Target 

% Vt = Target’s Speed 

Vt(l) = sqrt(Xt(2, 1 ) A 2+Xt(4, 1) A 2+Xt(6, 1 ) A 2); 

% Speed along the pitch and yaw LOS. Pitch and yaw closing speeds 
Vt_pitch(l) = Vt( 1 )*cos(LAMBD A(2, 1 )-GAMMA_t(2, 1 )); 

Vm_pitch(l) = Vm( 1 )*cos(LAMBDA(2, 1 )-GAMMA_m(2, 1 )); 

Vc_pitch(l) = Vm_pitch( 1 )*cos(GAMMA_m( 1 , 1 )-LAMB DA( 1,1)) 

- Vt_pitch( 1 )*cos(GAMMA_t( 1 . 1 )-LAMBD A( 1,1)); 

Vt_yaw(l) = Vt( l)*cos(GAMMA_t( 1,1)); 

Vm_yaw(l) = Vm(l)*cos(GAMMA_m(l,l)); 

Vc_yaw(l) = Vm_yaw(l)*cos(GAMMA_m(2,l)-LAMBDA(2,l)) 

- Vt_yaw( 1 )*cos(GAMMA_t(2, 1 )-LAMBD A(2, 1 )); 

Vc = [Vc_pitch(l) 0 

0 Vc_yaw(l)]; 

% SEEKER and GUIDANCE SYSTEM INITIAL CONDITIONS and INPUTS. 
% Seeker 
Xsk(:,l) = [0 
0 
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0 

0]; 

Usk(:,l) = LAMBDA(:,1); % Seeker input 
% Guidance System 
Xgs(:,l) = [0 
0 ]; 

Ugs(:,l) = N*Vc*[Xsk(2,l) % Guidance system input 
Xsk(4,l)]; 

% Initial conditions of the target’s acceleration 



a_tx(l) 


= 0 


a_ty(l) 


= 0 


a_tz(l) 


= 0 


a_t(l) 


= 0 


a_t_ pitch(l) 


= 0 


a_t_yaw(l) 


= 0; 



TETA_t(l) = 0; % angle between the acceleration vector and the yaw plane 
PHI_t(l) = 0; % yaw angle of the target’s acceleration vector 

% TIME 

% Time = Time vector 
TIME(l) = 0; 

% Tgo = Time to go 
Tgo(l) = R(4,l)/Vc_pitch(l); 

% SIMULATE THE SYSTEM 
for ti = 0:. 25:2 1.25 
for i = l:kmax-l 

% Calculate components of the missile’s pitch acceleration vector 
a_mx_pitch(i) = -(Xgs(l,i)*sin(LAMBDA(l,i))*cos(LAMBDA(2,i))); 
a_my_pitch(i) = -(Xgs(U)*sin(LAMBDA(ld))*sin(LAMBDA(2,i))); 
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a_mz_pitch(i) = Xgs(l,i)*cos(LAMBDA(l,i)); 

% Calculate components of the missile’s yaw acceleration vector 
a_mx_yaw(i) = -Xgs(2,i)*sin(LAMBDA(2,i)); 
a_my_yaw(i) = Xgs(2,i)*cos(LAMBDA(2,i)); 

% Compute overall missile acceleration components 
a_mx(i) = a_mx_pitch(i)+a_mx_yaw(i); 
a_my(i) - a_my_pitch(i)+a_my_yaw(i); 
a_mz(i) = a_mz_pitch(i); 
am(:,i) = [a_mx(i) 
a_my(i) 
a_mz(i)]; 

% target acceleration vector 
at(:,i) = [a_tx(i) 
a_ty(i) 
a_tz(i)]; 

% Compute magnitude of the missile’s acceleration 
a_m(i) - sqrt(a_mx(i) A 2 + a_my(i) A 2 + a_mz(i) A 2); 

% Generate target’s evasive maneuver (we assume that these accelerations, along 
% the three cartesian axis, are estimated using the missile’s image processing 
% capabilities) 

if TIME(i) >= ti/2 % target starts evasive maneuver 
a_tx(i+l) = 3*32.2*sin(GAMMA_t(2,i)); 
a_ty(i+l) = 4*32.2*cos(GAMMA_t(2,i)); 
a_tz(i+l) = 3*32.2*cos(GAMMA_t(l,i)); 
else 

a_tx(i+l) = 0.0; 
a_ty(i+l) = 0.0; 
a_tz(i+l) = 0.0; 
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end 



% Compute magnitude of the target’s acceleration 
a_t(i+l) = sqrt(a_tx(i+ 1 ) A 2 + a_ty(i+l) A 2 + a_tz(i+l) A 2); 

% Update missile states 

Xm(:,i+1) = PHIm*Xm(:,i)+DELm*am(:,i); 

% Update target states 

Xt(:,i+1) = PHIt*Xt(:,i)+DELt*at(:,i); 

% Update flight path angles 

GAMMA_m(:,i+l) = [atan2(Xm(6,i+l),sqrt(Xm(2,i+l) A 2+Xm(4,i+l) A 2)); 

atan2(Xm(4,i+l),Xm(2,i+l))]; 

GAMMA_t(:,i+l) = [atan2(Xt(6,i+l),sqrt(Xt(2,i+l) A 2+Xt(4,i+l) A 2)); 

atan2(Xt(4,i+ 1 ) ,Xt(2 ,i+ 1 ))] ; 

% Update seeker states 

Xsk(:,i+1) = PHIsk*Xsk(:,i)+DELsk*Usk(:,i); 

% Update Guidance System states 
Xgs(:,i+1) = PHIgs*Xgs(:,i)+DELgs*Ugs(:,i); 

% Limit yaw and pitch accelerations to 25 g’s 
if abs(Xgs(l,i+l)) > 805.0 
Xgs(l,i+1) = 805.0 *sign(Xgs(l,i+l)); 
end 

if abs(Xgs(2,i+l)) > 805.0 
Xgs(2,i+1) = 805.0 *sign(Xgs(2,i+l)); 
end 

% Update LOS angles 

LAMBD A(: ,i+ 1 ) = [atan2((Xt(5 ,i+ 1 )-Xm(5,i+ 1 )),sqrt((Xm( 1 ,i+ 1 )-Xt( 1 ,i+ 1 )) A 2 
+(Xm(3,i+l)-Xt(3,i+l)) A 2)); 

atan2((Xt(3,i+l)-Xm(3,i+l)),(abs(Xm(l,i+l)-Xt(l,i+l))))]; 
Usk(:,i+1) = LAMBDA(:,i+l); 
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LAMBDA_m(:,i+l) = [atan2(Xm(5,i+l),sqrt(Xm(l,i+l) A 2+Xm(3,i+l) A 2)); 

atan2(Xm(3,i+ 1 ),Xm( 1 ,i+ 1 ))] ; 

L AMBD A_t(: ,i+ 1 ) = [atan2(Xt(5 ,i+ 1 ),sqrt(Xt( 1 ,i+ 1 ) A 2+Xt(3,i+ 1 ) A 2)); 

atan2 (Xt(3 ,i+ 1 ),Xt( 1 , i+ 1 ))] ; 

% Update Range Information 

Rm(i+1) = sqrt(Xm(l,i+l) A 2 + Xm(3,i+1) A 2 + Xm(5,i+1) A 2); 

Rt(i+1) = sqrt(Xt(l,i+l) A 2 + Xt(3,i+1) A 2 + Xt(5,i+1) A 2); 

R(:,i+1) = [Xt(l,i+1)-Xm(l,i+1); 

Xt(3,i+1)-Xm(3,i+1); 

Xt(5,i+l)*Xm(5,i+l); 

sqrt((Xt( 1 ,i+ 1 )-Xm( 1 ,i+ 1 )) A 2+(Xt(3,i+ 1)-Xm(3,i+ 1 )) A 2 
+(Xt(5 ,i+ 1 )-Xm(5,i+ 1 )) A 2)]; 

% Update Velocity Information 

Vm(i+1) = sqrt(Xm(2,i+l) A 2+Xm(4,i+l) A 2+Xm(6,i+l) A 2); 

Vt(i+1) = sqrt(Xt(2,i+ 1 ) A 2+Xt(4,i+ 1) A 2+Xt(6,i+ 1) A 2); 

Vt_pitch(i+1) = Vt(i+l)*cos(LAMBDA(2,i+l)-GAMMA_t(2,i+l)); 
Vm_pitch(i+1) = Vm(i+l)*cos(LAMBDA(2,i+l)-GAMMA_m(2,i+l)); 
Vc_pitch(i+1) = Vm_pitch(i+l)*cos(GAMMA_m(l,i+l)-LAMBDA(l,i+l)) 
- Vt_pitch(i+ 1 )*cos(GAMMA_t( 1 ,i+ 1 )-LAMBD A( 1 ,i+ 1 )); 
Vt_yaw(i+1) = Vt(i+l)*cos(GAMMA_t(l,i+l)); 

Vm_yaw(i+1) = Vm(i+l)*cos(GAMMA_m(l,i+l)); 

Vc_yaw(i+1) = Vm_yaw(i+l)*cos(GAMMA_m(2,i+l)-LAMBDA(2,i+l)) 
-Vt_yaw(i+l)*cos(GAMMA_t(2,i+l)-LAMBDA(2,i+l)); 
Vc = [Vc_pitch(i+1) 0 

0 Vc_yaw(i+1)]; 

% Calculate angles of the target’s acceleration 
TETA_t(i+l) = atan2(a_tz(i+l),sqrt(a_tx(i+l) A 2+a_ty(i+l) A 2)); 

PHI_t(i+l) = atan2(a_ty(i+l),a_tx(i+l)); 
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% Calculate the components of the target’s acceleration normal to the LOS 
a_t_pitch(i+l) = -a_t(i+l)*cos(LAMBDA(2,i+l)-PHI_t(i+l)) 

*sin(LAMBDA( 1 ,i+ l)-TETA_t(i+ 1 )); 

a_t_yaw(i+l) = -a_t(i+l)*cos(TETA_t(i+l))*sin(LAMBDA(2,i+l)-PHI_t(i+l)); 
% Update guidance system input 
Ugs(:,i+1) = N*(Vc*[Xsk(2,i+l);Xsk(4,i+l)] 

+.5*[a_t_pitch(i+l);a_t_yaw(i+l)]); 

% Update Time/time to go 
TIME(i+l) =TIME(i)+dt; 

Tgo(i+l) = R(4,i+1)/Vc_pitch(i+1); 

% Check for closest point 
if (R(4,i) < R(4,i+l)),break,end 
end; 

% Save information for plotting and evaluation 
Rl(4*ti+1) = R(4,i); % miss distance 
Ti(4*ti+1) = ti/2;% starting time of evasive maneuver (EM) 
tgo(4*ti+l) = i*dt-ti/2; % time to go until end of flight 
if ti == 12.0 % Record information for a target that 

% initialized the evasive maneuver 6 sec 
% after the missile entered into the terminal 
TGO = tgo(49); % phase of flight 

Xseeker = Xsk(:,l:i); 

Xgsys = Xgs(:,l:i); 
lambda_m = LAMBDA_m(:,l:i); 
lambda_t = LAMBDA_t(:,l:i); 
lambda = LAMBDA(:,l:i); 
gamma_m = GAMMA_m(:,l:i); 
gamma_t = GAMMA_t(:,l:i); 
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r = R(:,l:i); 
vm = Vm(l:i); 
vt = Vt(l:i); 

vm_pitch = Vm_pitch(l:i); 
vt_pitch = V t_pitch( 1 : i); 
vm_yaw = Vm_yaw(l:i); 
vt_yaw = Vt_yaw(l:i); 
vc_pitch = Vc_pitch(l:i); 
vc_yaw = Vc_yaw(l:i); 
tGO = Tgo(l:i); 
a_M = am(:,l:i); 
a_T = at(:,l:i); 

A_t = a_t(l:i); 

A_m = a_m(l:i); 
time = TIME(l:i); 
end 

clear R; 

R(:,l) = [Xt(l,l)-Xm(l,l); 

Xt(3,l)-Xm(3,l); 

Xt(5,l)-Xm(5,l); 

sqrt((Xt(l,l)-Xm(l,l)) A 2+(Xt(3,l)-Xm(3,l)) A 2+(Xt(5,l)-Xm(5,l)) A 2)]; 

end; 

save thesis2a343 R1 tgo Ti tGO missile target TGO Xseeker Xgsys lambda_m 
lambda_t lambda gamma_m gamma_t r vm vt vm_pitch vm_yaw vt_pitch 
vt_yaw vc_pitch vc_yaw a_M a_T time A_t A_m 
PLOTS 

% Miss distance information 

plot(Ti,R 1 ) ,title( ‘MISS DISTANCE vs INITIAL TIME, APROPNAV’) 
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xlabel(‘ INITIAL TIME - SEC’),ylabel(‘MISS - FEET’) 

print -dps Rlaal 

Ipstoepsi Rlaal.ps Rlaal.epsi 

pause, clg 

plot(tgo,Rl),title(‘MISS DISTANCE vs TIME TO GO, APROPNAV’) 

xlabel(‘TIME TO GO - SEC’),ylabel(‘MISS - FEET’) 

print -dps Rlbal 

Ipstoepsi Rlbal.ps Rlbal.epsi 

pause, clg 

% Missile acceleration information 

plot(time, A_m) ,title( ‘ MISS ILE ACCELERATION MAGNITUDE vs TIME, 

APROPNAV’) 

xlabel(‘TTME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps A_mal 
Ipstoepsi A_mal.ps A_mal.epsi 
pause.clg 

plot(time,Xgsys(l,:)),title(‘MISSILE PITCH ACCELERATION vs TIME, 

APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps Xgsyslal 
Ipstoepsi Xgsyslal.ps Xgsyslal. epsi 
pause.clg 

plot(time,Xgsys(2,:)),title(‘MISSILE YAW ACCELERATION vs TIME, 

APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps Xgsys2al 
Ipstoepsi Xgsys2al.ps Xgsys2al.epsi 
pause.clg 
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% Target acceleration information 

plot(time,A_t),title(‘TARGET ACCELERATION MAGNITUDE vs TIME, 
APROPNAV’) 

xlabel(‘TIME - S EC ’ ) , y label( ‘ FEET/S EC A 2 ’ ) 

print -dps A_tal 

Ipstoepsi A_tal.ps A_tal.epsi 

pause,clg 

% Seeker pitch and yaw angles 

plot(time,Xseeker(l,:)),title(‘SEEKEK PITCH ANGLE vs TIME, APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD’) 

print -dps Xseekerlal 

Ipstoepsi Xseekerlal.ps Xseekerlal. epsi 

pause,clg 

plot(time,Xseeker(3,:)),title(‘SEEKER YAW ANGLE vs TIME, APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD’) 

print -dps Xseeker2al 

Ipstoepsi Xseeker2al.ps Xseeker2al.epsi 

pause,clg 

plot(time,Xseeker(2,:)),title(‘SEEKER PITCH ANGLE RATE vs TIME, 

APROPNAV’) 

xlabel(‘TIME - S EC ’),ylabel(‘ RAD/SEC’) 
print -dps Xseeker3al 
Ipstoepsi Xseeker3al.ps Xseeker3al.epsi 
pause.clg 

plot(time,Xseeker(4,:)),title(‘SEEKER YAW ANGLE RATE vs TIME, 

APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD/SEC’) 
print -dps Xseeker4al 
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Ipstoepsi Xseeker4al.ps Xseeker4al.epsi 
pause,clg 

% Range information 

plot(time,r(4,:)),title(‘RANGE vs TIME, APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET’) 

print -dps ral 

Ipstoepsi ral.ps ral.easi 

pause,clg 

% Missile velocity information 

plot(time,vm),title(‘MISSILE SPEED vs TIME, APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vmal 

Ipstoepsi vmal.ps vmal.epsi 

pause,clg 

% Target velocity information 

plot(time,vt),title(‘TARGET SPEED vs TIME, APROPNAV’) 

xlabel(‘TTME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vtal 

Ipstoepsi vtal.ps vtal.epsi 

pause,clg 

% Closing velocity information 

plot( time, vc_pitch) , title( ‘ PIT C H CLOSING SPEED, APROPNAV’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vclal 

Ipstoepsi vclal.ps vclal.epsi 

pause, clg 

plot( time,vc_yaw),title( ‘ Y AW CLOSING SPEED vs TIME, APROPNAV’) 
xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 
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print -dps vc2al 

Ipstoepsi vc2al.ps vc2al.epsi 

pause, clg 
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APPENDIX C - MISSILE/TARGET THREE DIMENSIONAL 
SIMULATION USING OPTIMAL GUIDANCE 



% Written by: Rui Manuel Alves Francisco 
% Date: 09 December 1992 

% This Program simulates the terminal phase of a 3-D missile/target 
% engagement using optimal guidance. 

clear 

clg 

% DEFINE CONSTANTS 
dt = .01; % Sampling time 
Tf = 100; % maximum simulation time 
kmax = Tf/dt+1; 

% DEFINE STATE EQUATIONS 
% Missile 

% Xm = [xm = missile's x coordinate 
% xmd = missile's speed (x coordinate) 

% ym = missile's y coordinate 

% ymd = missile's speed (y coordinate) 

% zm = missile's z coordinate 

% zmd = missile's speed (z coordinate)] 

Am = [0 1 0 0 0 0 
000000 
000100 
000000 
000001 
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00 00 00]; 

Bm = [000 
100 
000 
0 10 
000 
0 0 1 ]; 

% Target 

% Xt = [xt = target’s x coordinate 
% xtd = target’s speed (x coordinate) 
% yt = target’s y coordinate 

% ytd = target’s speed (y coordinate) 

% zt = target’s z coordinate 

% ztd = target’s speed (z coordinate)] 

At = [0 1 0 0 0 0 
000000 
000100 
000000 
000001 
0 0 00 0 0 ]; 

Bt = [0 0 0 
100 
000 
0 1 0 
000 
0 0 1]; 

% Seeker (Radar) 
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% Xsk = [beta_pitch = seeker’s pitch angle 
% betad_pitch = seeker’s pitch angle rate 

% beta_yaw = seeker’s yaw angle 

% betad_yaw = seeker’s yaw angle rate] 



Ask = [ 0 1 0 0 

-100 -20 0 0 

0 0 0 1 

0 0 -100 - 20 ]; 

Bsk = [ 0 0 

100 0 
0 0 
0 100 ]; 

% Guidance System 

% Xgs = [a_m_pitch = missile’s pitch acceleration 
% a_m_yaw = missile’s yaw acceleration] 

Ags = [-10 

0-1]; 

Bgs = [10 

0 1 ]; 

% INITIALIZE STATE VARIABLES (when the missile enters into the terminal 
phase of flight) 



% Missile 



Xm(:,l) = [ 0 

2828 



% The missile is in a collision triangle 
% with the target when the missile enters into 
% the terminal phase of flight 



0 



1000 



0 
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47.1339]; 



% Target 
Xt(:,l) = [30000 
0 
0 

1000 

500 

0 ]; 

% DISCRETE REPRESENTATION 
[PHIm,DELm] = c2d(Am,Bm,dt); 

[PHIt,DELt] = c2d(At,Bt,dt); 

[PHIsk.DELsk] = c2d(Ask,Bsk,dt); 

[PHIgs,DELgs] = c2d(Ags,Bgs,dt); 

% LINE OF SIGHT (LOS) INFORMATION. INITIAL CONDITIONS. 
% Missile 

% LAMBDA_m = Missile’s LOS from the global coordinate system 
% LAMBDA_m = [LAMBDA_m_pitch = Missile’s pitch LOS angle 
% LAMBDA_m_yaw = Missile’s yaw LOS angle] 

LAMBDA_m(:,l) = [atan2(Xm(5,l),sqrt(Xm(Ll) A 2+Xm(3,l) A 2)); 

atan2(Xm(3 , 1 ),Xm( 1,1))]; 

% Target 

% LAMBDA_t = Target’s LOS from the global coordinate system 
% LAMBDA_t = [LAMB DA_t_p itch = Target’s pitch LOS angle 
% LAMB D A_t_y aw = Target’s yaw LOS angle] 

LAMBDA_t(:,l) = [atan2(Xt(5, 1 ),sqrt(Xt( 1 , 1 ) A 2+Xt(3, 1) A 2)); 

atan2(Xt(3,l) ,Xt( 1,1))]; 

% LOS from Missile to Target 
% LAMBDA = LOS from Missile to Target. 
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% LAMBDA = [LAMB DA_p itch = LOS angle in pitch 
% LAMBDA_yaw = LOS angle in yaw]; 

LAMBDA(:,1) = [atan2((Xt(5,l)-Xm(5,l)),sqrt((Xt(l,l)-Xm(l,l)) A 2 
+(Xt(3,l)-Xm(3,l)) A 2)); 

atan2((Xt(3, 1 )-Xm(3, 1 )),abs(Xt( 1 , 1 )-Xm( 1,1)))]; 

% MISSILE and TARGET FLIGHT PATH ANGLES INFORMATION 
% Missile 

%GAMMA_m = [GAMMA_m_pitch = Missile’s flight path angle in pitch 
% GAMMA_m_yaw = Missile’s flight path angle in yaw] 

GAMMA_ra(:,l) = [atan2(Xm(6,l),sqrt(Xm(2,l) A 2+Xm(4,l) A 2)); 
atan2(Xm(4,l),Xm(2,l))]; 

% Target 

%GAMMA_t = [GAMMA_t_pitch = Target’s flight path angle in pitch 
% GAMMA_t_YAW = Target’s flight path angle in yaw] 

GAMMA_t(:,l) = [atan2(Xt(6,l),sqrt(Xt(2,l) A 2+Xt(4,l) A 2)); 

atan2(Xt(4, 1 ) ,Xt(2, 1 ))] ; 

% RANGE INFORMATION 
% Missile 

% Rm = Missile’s range 

Rm(l) = sqrt(Xm(l,l) A 2 + Xm(3,l) A 2 + Xm(5,l) A 2); 

% Target 

% Rt = Target’s range 

Rt(l) = sqrt(Xt(l,l) A 2 + Xt(3,l) A 2 + Xt(5,l) A 2); 

% Missile/Target relative distance 
% R = [Rmtx = Missile/Target x coordinate range 
% Rmty = Missile/Target y coordinate range 

% Rmtz = Missile/Target z coordinate range 

% Rmt = Missile/Target relative distance(Rmt=sqrt(Rmtx A 2+Rmty A 2+Rmtz A 2))] 
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R(:,l) = [Xt(l,l)-Xm(l,l) 

Xt(3,l)-Xm(3.1) 

Xt(5,l)-Xm(5,l) 

sqrt((Xt(l,l)-Xm(l,l)) A 2+(Xt(3,l)-Xm(3,l)) A 2+(Xt(5,l)-Xm(5,l)) A 2)]; 

% VELOCITY INFORMATION 
% Missile 

% Vm = Missile’s Speed 

Vm(l) = sqrt(Xm(2,l) A 2+Xm(4,l) A 2+Xm(6,l) A 2); 

% Target 

% Vt = Target’s Speed 

Vt(l) = sqrt(Xt(2 , 1 ) A 2+Xt(4, 1 ) A 2+Xt(6, 1 ) A 2); 

% Speed along the pitch and yaw LOS. Pitch and yaw closing speeds 
Vt_pitch(l) = Vt(l)*cos(LAMBDA(2,l)-GAMMA_t(2,l)); 

Vm_pitch(l) = Vm(l)*cos(LAMBDA(2,l)-GAMMA_m(2,l)); 

Vc_pitch(l) = Vm_pitch( 1 )*cos(G AMMA_m( 1 , 1 )-LAMB D A( 1,1)) 
-Vt_pitch(l)*cos(GAMMA_t(l,l)-LAMBDA(l,l)); 
Vt_yaw(l) = Vt(l)*cos(GAMMA_t(l,l)); 

Vm_yaw(l) = Vm(l)*cos(GAMMA_m(l,l)); 

Vc_yaw(l) = Vm_yaw(l)*cos(GAMMA_m(2,l)-LAMBDA(2,l)) 
*Vt_yaw(l)*cos(GAMMA_t(2,I)-LAMBDA(2,l)); 

Vc = [Vc_pitch(l) 0 

0 Vc_yaw(l)]; 

% Tgo = Time to go 
Tgo(l) = R(4,l)/Vc_pitch(l); 

% Optimal guidance coefficients 
k= Tgo(l); 

n = (6*k A 2*(exp(-k)-l+k))/(2*k A 3+3+6*k-6*k A 2-12*k*exp(-k)-3*exp(-2*k)); 
N = [n 0 
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0 n]; 

% SEEKER and GUIDANCE SYSTEM INITIAL CONDITIONS and INPUTS. 
% Seeker 
Xsk(:,l) = [0 
0 
0 

0 ]; 

Usk(:,l) = LAMBDA(:,1); % Seeker input 
% Guidance System 
Xgs(:,l) = [0 
0 ]; 

Ugs(:,l) = N*Vc*[Xsk(2,l) % Guidance system input 
Xsk(4,l)]; 

% Initial conditions of the target’s acceleration 



a_tx(l) 


= 0; 


a_ty(l) 


= 0; 


a_tz(l) 


= 0; 


a_t(l) 


= 0; 


a _ t_pitch(l) 


= 0; 


a_t_yaw(l) 


= 0; 


TETA_t(l) 


= 0; % angle between the acceleration vector and the 


PHI_t(l) 


= 0; % yaw angle of the target’s acceleration vector 


% Initial conditions of the missile’s acceleration 


a_m_pitch(l) 


= 0; 


a_m_yaw(l) 


= 0; 



% TIME 

% Time = Time vector 
TTME(l) = 0; 
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% SIMULATE THE SYSTEM 
for ti = 0:. 25:2 1.25 
for i = l:kmax-l 

% Calculate components of the missile’s pitch acceleration vector 
a_mx_pitch(i) = -(Xgs(l,i)*sin(LAMBDA(l,i))*cos(LAMBDA(2,i))); 
a_my_pitch(i) = -(Xgs(l,i)*sin(LAMBDA(l,i))*sin(LAMBDA(2,i))); 
a_mz_pitch(i) = Xgs(l,i)*cos(LAMBDA(l,i)); 

% Calculate missile’s yaw acceleration vector components 
a_mx_yaw(i) = -Xgs(2,i)*sin(LAMBDA(2,i)); 
a_my_yaw(i) = Xgs(2,i)*cos(LAMBDA(2,i)); 

% Compute overall missile acceleration 
a_mx(i) = a_mx_pitch(i)+a_mx_yaw(i); 
a_my(i) = a_my_pitch(i)+a_my_yaw(i); 
a_mz(i) = a_mz_pitch(i); 
am(:,i) = [a_mx(i) 
a_my(i) 
a_mz(i)]; 

% target acceleration vector 
at(:,i) = [a_tx(i) 
a_ty(i) 
a_tz(i)]; 

% Compute magnitude of the missile’s acceleration 
a_m(i) = sqrt(a_mx(i) A 2 + a_my(i) A 2 + a_mz(i) A 2); 

% Generate target’s evasive maneuver (we assume that these accelerations, along 
% the three cartesian axis, are estimated using the missile’s image processing 
% capabilities) 

if TIME(i) >= ti/2 % target starts evasive maneuver 
a_tx(i+l) = 3*32.2*sin(GAMMA_t(2,i)); 
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a_ty(i+l) = 4*32.2*cos(GAMMA_t(2,i)); 
a_tz(i+l) = 3*32.2*cos(GAMMA_t(l,i)); 
else 

a_tx(i+l) = 0.0; 
a_ty(i+l) = 0.0; 
a_tz(i+l) = 0.0; 
end 

% Compute magnitude of the target’s acceleration 
a_t(i+l) = sqrt(a_tx(i+l) A 2 + a_ty(i+l) A 2 + a_tz(i+l) A 2); 

% Update missile states 

Xm(:,i+1) = PHIm*Xm(:,i)+DELm*am(:,i); 

% Update target states 

Xt(:,i+1) = PHIt*Xt(:,i)+DELt*at(:,i); 

% Update flight path angles 

GAMMA_m(:,i+l) = [atan2(Xm(6,i+l),sqrt(Xm(2,i+l) A 2+Xm(4,i+l) A 2)); 

atan2(Xm(4,i+ 1 ),Xm(2,i+ 1 ))] ; 

GAMMA_t(:,i+l) = [atan2(Xt(6,i+l),sqrt(Xt(2,i+l) A 2+Xt(4,i+l) A 2)); 

atan2 (Xt(4 ,i+ 1 ) ,Xt(2 ,i+ 1 ))] ; 

% Update seeker states 

Xsk(:,i+1) = Pfflsk*Xsk(:,i)+DELsk*Usk(:,i); 

% Update Guidance System states 
Xgs(:,i+1) = Pfflgs*Xgs(:,i)+DELgs*Ugs(:,i); 

% Limit yaw and pitch accelerations to 25 g’s 
if abs(Xgs(l,i+l))> 805.0 
Xgs(l,i+1) = 805.0 *sign(Xgs(l,i+l)); 
end 

if abs(Xgs(2,i+l)) > 805.0 
Xgs(2,i+1) = 805.0 *sign(Xgs(2,i+l)); 
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end 



% Update LOS angles 

LAMBDA(:,i+ 1 ) = [atan2((Xt(5,i+ 1 )-Xm(5 ,i+ 1 )),sqrt((Xm( 1 ,i+ 1 )-Xt( 1 ,i+ 1)) A 2 
+(Xm(3,i+l)-Xt(3,i+l)) A 2)); 

atan2((Xt(3,i+l)-Xm(3,i+l)),(abs(Xm(l,i+l)-Xt(l,i+l ))))]; 
Usk(:,i+1) = LAMBDA(:,i+l); 

LAMBDA_m(:,i+l) = [atan2(Xm(5,i+l),sqrt(Xm(l,i+l) A 2+Xm(3,i+l) A 2)); 

atan2(Xm(3,i+l),Xm(l,i+l))]; 

LAMB D A_t(: 4 + 1) = [atan2(Xt(5 ,i+ 1 ),sqrt(Xt( 1 ,i+ 1 ) A 2+Xt(3 ,i+ 1 ) A 2)); 

atan2(Xt(3 ,i+ 1 ),Xt( 1 ,i+ 1 ))] ; 

% Update Range Information 

Rm(i+1) = sqrt(Xm(l,i+l) A 2 + Xm(3,i+1) A 2 + Xm(5,i+1) A 2); 

Rt(i+1) = sqrt(Xt(l,i+l) A 2 + Xt(3,i+1) A 2 + Xt(5,i+1) A 2); 

R(:,i+1) = [Xt( 1 ,i+l )-Xm( 1 ,i+ 1 ); 

Xt(3,i+ 1 )-Xm(3,i+ 1 ); 

Xt(5 ,i+ 1 )-Xm(5 ,i+ 1 ); 

sqrt((Xt(l,i+l)-Xm(l,i+l)) A 2+(Xt(3,i+l)-Xm(3,i+l)) A 2 

+(Xt(5,i+l)-Xm(5,i+l)) A 2)]; 

% Update Velocity Information 

Vm(i+1) = sqrt(Xm(2,i+ 1) A 2+Xm(4,i+ 1) A 2+Xm(6,i+ 1 ) A 2); 

Vt(i+1) = sqrt(Xt(2,i+ 1 ) A 2+Xt(4,i+ 1) A 2+Xt(6,i+ 1) A 2); 

Vt_pitch(i+1) = Vt(i+l)*cos(LAMBDA(2,i+l)-GAMMA_t(2,i+l)); 
Vm_pitch(i+1) = Vm(i+l)*cos(LAMBDA(2,i+l)-GAMMA_m(2,i+l)); 
Vc_pitch(i+1) = Vm_pitch(i+1 )*cos(GAMMA_m( 1 ,i+ 1 )-LAMBD A( 1 ,i+ 1 )) 
-Vt_pitch(i+l)*cos(GAMMA_t(l,i+l)-LAMBDA(l,i+l)); 
Vt_yaw(i+1) = Vt(i+l)*cos(GAMMA_t(l,i+l)); 

Vm_yaw(i+1) = Vm(i+l)*cos(GAMMA_m(l,i+l)); 

Vc_yaw(i+1) = Vm_yaw(i+l)*cos(GAMMA_m(2,i+l)-LAMBDA(2,i+l)) 
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-Vt_yaw(i+l)*cos(GAMMA_t(2,i+l)-LAMBDA(2,i+l)); 

Vc = [Vc_pitch(i+1) 0 

0 Vc_yaw(i+1)]; 

% Time to go 

Tgo(i+l) = R(4,i-t- 1 )/Vc_pitch(i-»-l ); 

% Calculate angles of the target’s acceleration 
TETA_t(i+l) = atan2(a_tz(i+l),sqrt(a_tx(i+l) A 2+a_ty(i+l) A 2)); 

PHI_t(i+l) = atan2(a_ty(i+l),a_tx(i+l)); 

% Calculate the components of the target’s acceleration normal to the LOS 
a_t_pitch(i+l) = -a_t(i+l)*cos(LAMBDA(2,i+l)-PHI_t(i+l)) 

*sin(LAMB D A( 1 ,i+ 1 )-TETA_t(i+l )); 

a_t_yaw(i+l) = -a_t(i+l)*cos(TETA_t(i+l))*sin(LAMBDA(2,i+l)-PHI_t(i+l)); 
% Update optimal guidance coefficients 
k = Tgo(i+l); 

n = (6*k A 2*(exp(-k)-l+k))/(2*k A 3+3+6*k-6*k A 2-l2*k*exp(-k)-3*exp(-2*k)); 
N=[n 0 
0 n]; 

% Components of the Missile’s acceleration normal to the pitch and yaw LOS 
a_m_pitch(i+l) = Xgs(l,i+1); 
a_m_yaw(i+l) = Xgs(2,i+1); 

% Update guidance system input 
Ugs(:,i+1) = N*(Vc*[Xsk(2,i+l);Xsk(4,i+l)] 

+.5*[a_t_pitch(i+l);a_t_yaw(i+l)]) 

-(l/k A 2)*(exp(-k)+k-l)*[a_m_pitch(i+l);a_m_yaw(i+l)]); 

% Update Time 
TIME(i+l) = TIME(i)+dt; 

% Check for closest point 
if (R(4,i) < R(4,i+l)),break,end 
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end; 

% Save information for plotting and evaluation 

Rl(4*ti+1) = R(4,i); % miss distance 

Ti(4*ti+1) = ti/2;% starting time of evasive maneuver (EM) 

tgo(4*ti+l) = i*dt-ti/2; % time to go until end of flight 

if ti == 12.0 % Record information for a target that 

% initialized the evasive maneuver 6 sec after 
% the missile entered into the terminal phase 
TGO = tgo(49); % of flight 

Xseeker = Xsk(:,l:i); 

Xgsys = Xgs(:,l:i); 
lambda_m = LAMBDA_m(:,l:i); 
lambda_t = LAMBDA_t(:,l:i); 
lambda = LAMBDA(:,l:i); 
gamma_m = GAMMA_m(:,l;i); 
gamma_t = GAMMA_t(:,l;i); 
r = R(:,l:i); 
vm = Vm(l:i); 
vt = Vt(l;i); 

vm_pitch = Vm_pitch(l:i); 
vt_pitch = Vt_pitch(l;i); 
vm_yaw = Vm_yaw(l:i); 
vt_yaw = Vt_yaw(l:i); 
vc_pitch = Vc_pitch(l:i); 
vc_yaw = Vc_yaw(l:i); 
tGO =Tgo(l:i); 
a_M = am(:,l:i); 
a_T = at(:,l:i); 



143 



A_t = a_t(l:i); 

A_m = a_m(l:i); 
time = TIME(l:i); 
end 

clear R; 

R(:,l) = [Xt(l,l)-Xm(l,l); 

Xt(3,l)-Xm(3,l); 

Xt(5,l)-Xm(5,l); 

sqrt((Xt(l,l)-Xm(l,l)) A 2+(Xt(3,l)-Xm(3,l)) A 2+(Xt(5,l)-Xm(5,l)) A 2)]; 

end; 

save thesis3o343 R1 tgo Ti tGO missile target TGO Xseeker Xgsys lambda_m 
lambda_t lambda gamma_m gamma_t r vm vt vm_pitch vm_yaw vt_pitch 
vt_yaw vc_pitch vc_yaw a_M a_T time A_t A_m 
PLOTS 

% Miss distance information 

plot(Ti,R 1 ),tide( ‘MISS DISTANCE vs. INITIAL TIME, OPTIMAL’) 

xlabel( ‘INITIAL TIME - SEC’),ylabel(‘MISS - FEET’) 

print -dps Rlaol 

Ipstoepsi Rlaol.ps Rlaol.epsi 

pause.clg 

plot(tgo,Rl),title(‘MISS DISTANCE vs TIME TO GO, OPTIMAL’) 

xlabel(‘TIME TO GO - SEC’),ylabel(‘MISS - FEET’) 

print -dps Rlbol 

Ipstoepsi Rlbol.ps Rlbol.epsi 

pause, clg 

% Missile acceleration information 

plot(time,A_m),title(‘MISSILE ACCELERATION MAGNITUDE vs TIME, 

OPTIMAL’) 
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xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps A_mol 
Ipstoepsi A_mol.ps A_mol.epsi 
pause.clg 

plot(time,Xgsys(l,:)),title(‘MISSILE PITCH ACCELERATION vs TIME, 

OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps Xgsyslol 
Ipstoepsi Xgsyslol.ps Xgsyslol.epsi 
pause,clg 

plot(time,Xgsys(2,:)),title( ‘MISSILE YAW ACCELERATION vs TIME, 

OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 
print -dps Xgsys2ol 
Ipstoepsi Xgsys2ol.ps Xgsys2ol.epsi 
pause, clg 

% Target acceleration information 

plot(time,A_t),title(‘TARGET ACCELERATION MAGNITUDE vs TIME, 
OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC A 2’) 

print -dps A_tol 

Ipstoepsi A_tol.ps A_tol.epsi 

pause, clg 

% Seeker pitch and yaw angles 

plot(time,Xseeker(l,:)),title(‘SEEKER PITCH ANGLE vs TIME, OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD’) 

print -dps Xseekerlol 

Ipstoepsi Xseekerlol.ps Xseekerlol.epsi 
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pause, clg 

plot(time,Xseeker(3,:)),title(‘SEEKER YAW ANGLE vs TIME, OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD’) 

print -dps Xseeker2ol 

Ipstoepsi Xseeker2ol.ps Xseeker2ol.epsi 

pause, clg 

plot(time,Xseeker(2,:)),title(‘SEEKER PITCH ANGLE RATE vs TIME, 

OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD/SEC’) 
print -dps Xseeker3ol 
Ipstoepsi Xseeker3ol.ps Xseeker3ol.epsi 
pause,clg 

plot(time,Xseeker(4,:)),title(‘SEEKER YAW ANGLE RATE vs TIME, 

OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘RAD/SEC’) 
print -dps Xseeker4ol 
Ipstoepsi Xseeker4ol.ps Xseeker4ol.epsi 
pause,clg 

% Range information 

plot(time/(4,:)),title(‘RANGE vs TIME, OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET’) 

print -dps rol 

Ipstoepsi rol.ps rol.easi 

pause,clg 

% Missile velocity information 

plot(time,vm),title(‘MISSILE SPEED vs TIME, OPTIMAL’) 
xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 
print -dps vmol 
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Ipstoepsi vmol.ps vmol.epsi 
pause, clg 

% Target velocity information 

plot(time,vt),title( ‘TARGET SPEED vs TIME, OPTIMAL’) 

xlabel(‘TEME - SEC’),ylabel(‘FEET/SEC') 

print -dps vtol 

Ipstoepsi vtol.ps vtol.epsi 

pause,clg 

% Closing velocity information 

plot(time,vc_pitch),title(‘PITCH CLOSING SPEED, OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vclol 

Ipstoepsi vclol.ps vclol. epsi 

pause, clg 

plot(time,vc_yaw),title(‘YAW CLOSING SPEED vs TIME, OPTIMAL’) 

xlabel(‘TIME - SEC’),ylabel(‘FEET/SEC’) 

print -dps vc2ol 

Ipstoepsi vc2ol.ps vc2ol.epsi 

pause, clg 
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